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Abstract 

When  using  a  thin/film  gage  to  measure  surface  heat  flux,  one  typically  re¬ 
duces  the  data  for  surface  temperature  to  surface  heat  flux  using  a  series  solution  for 
one-dimensional  heat  transfer  in  an  initially  isothermal,  semi-infinite  solid.  How¬ 
ever,  the  gage  may  not  behave  as  an  initially  isothermal,  semi-infinite  solid  due 
to  multidimensional  heat  transfer  and  electrical  preheating  of  the  gage  when  the 
instrumentation  is  tinned  on. 


To  evaluate  the  accuracy  of  the  series  solution  for  use  with  thin-film  gages, 
the  heat  transfer  in  a  gage  was  numerically  simulated  using  a  two-dimensional, 
finite-difference  model.  The  actual  geometry  of  the  probe  was  simplified  to  reduce 
the  heat  transfer  to  two  dimensions.  The  simulation  produced  surface  temperatures 
which  were  used  in  the  series  solution  to  find  estimates  of  surface  heat  flux.  The 
heat  fluxes  from  the  simulation  and  the  series  solution  were  then  compared  to 
evaluate  the  accuracy  of  the  series  solution. 

The  analysis  provides  good  insight  into  the  causes  of  inaccuracies  when  using 
the  series  solution.  It  also  provides  some  quantitative  results  which  may  be  helpful 
for  estimating  errors  in  actual  laboratory  use.  /K  ,  ,  X 
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A  NUMERICAL  INVESTIGATION  OF  THIN-FILM 
HEAT  TRANSFER  GAGES 


I.  Introduction 

1.1  Background  and  Motivation 

A  common  experimental  technique  for  measuring  surface  heat  flux  in  the 
laboratory  employs  a  thin-film  gage.  A  thin-film  gage  is  essentially  a  small  quartz 
cylinder  which  has  a  thin,  narrow  strip  of  platinum — the  thin-film  sensor — plated 
on  its  front  face.  Two  small  wires  connected  to  opposite  ends  of  the  thin-film 
sensor  run  along  the  side  of  the  quartz  cylinder  to  connect  the  sensor  to  recording 
instrumentation.  Figure  1  is  a  diagram  of  a  thin-film  gage. 

The  thin-film  gage  does  not  directly  measure  surface  heat  flux.  Instead, 
it  measures  the  change  in  temperature  of  the  surface  of  the  quartz  cylinder  at 
the  location  of  the  thin-film  sensor  as  a  fuction  of  time.  Because  the  platinum 
film  is  thin  and  has  a  relatively  high  thermal  conductivity,  the  temperature  of 
the  thin-film  sensor  should  be  the  same  as  the  temperature  of  the  quartz  cylinder 
directly  beneath  it.  The  resistance  of  the  thin-film  sensor  is  a  known  function 
of  temperature.  Each  gage  is  calibrated  to  determine  the  change  in  resistance 
as  a  function  of  change  in  temperature  above  some  reference  temperature  [5:p.9]. 
As  the  surface  temperature  of  the  cylinder  changes  in  response  to  a  disturbance, 
the  electrical  resistance  of  the  thin-film  sensor  will  change.  Then,  if  the  gage 
is  connected  in  one  leg  of  a  properly  balanced  Wheatstone  bridge,  the  output 
voltage  of  the  Wheatstone  bridge  will  be  proportional  to  the  change  in  resistance 
of  the  thin-film  sensor.  Thus,  the  output  voltage  of  the  Wheatstone  bridge  circuit 
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indirectly  measures  the  change  in  surface  temperature  of  the  area  of  the  quartz 
cylinder  directly  beneath  the  thin-film  sensor. 

Balancing  the  Wheatstone  bridge  requires  consideration  of  the  electrical  heat¬ 
ing  in  the  thin-film  sensor.  When  the  input  voltage  is  applied  to  the  bridge,  a  small 
current  will  flow  through  each  leg  of  the  bridge  circuit  including  the  thin-film  sen¬ 
sor.  This  small  current  will  cause  some  electrical  heating  in  the  thin-film  sensor 
changing  the  resistance  of  the  sensor  and  unbalancing  the  bridge.  Therefore,  the 
circuit  must  be  turned  on  and  allowed  to  warm  up  until  the  transients  settle  down 
before  the  bridge  can  be  balanced. 

To  measure  the  surface  heat  flux  at  the  front  surface  of  a  test  specimen,  a 
thin-film  gage  is  embedded  in  the  specimen.  The  gage  is  held  in  place  with  a  filler 
material  which  provides  a  sealed  fit.  The  front  surface  of  the  gage  must  be  flush 
with  the  front  surface  of  the  test  specimen. 

To  get  results  for  surface  heat  flux  using  the  thin-film  gage,  one  must  convert 
the  data  for  change  in  surface  temperature  as  a  function  of  time  into  values  for 
surface  heat  flux  as  a  function  of  time.  In  practice,  the  data  is  reduced  using  a 
series  solution  which  assumes  heat  transfer  in  an  initially  isothermal,  semi-infinite 
solid.  In  many  ways,  however,  the  gage  does  not  behave  as  an  initially  isothermal, 
semi-infinite  solid  during  the  time  when  surface  temperature  data  is  being  collected. 
First  of  all,  the  electrical  heating  in  the  thin-film  sensor  that  occurs  while  the  bridge 
is  being  balanced  produces  a  non-uniform,  initial  temperature  distribution  in  the 
quartz  cylinder.  Secondly,  the  electrical  heating  at  any  time  generates  non-uniform 
heat  flow  into  the  quartz  cylinder  which  induces  a  three-dimensional  temperature 
profile  and  heat  flow  pattern.  Thirdly,  heat  flows  radially  in  the  presence  of  the 
disturbance  because  the  gage  is  not  thermally  isolated  from  its  surroundings  and 
the  thermal  properties  of  the  surrounding  materials  are  not  necessarily  identical 
to  those  of  the  quartz  cylinder.  For  the  results  for  external  surface  heat  flux  to 
be  accurate,  the  effect  of  the  multidimensional  heat  transfer  on  the  series  solution 


must  be  small. 

Ideally,  one  hopes  to  measure  the  external  surface  heat  flux  that  would  be 
present  in  the  test  specimen  if  the  gage  were  not  there.  The  external  surface 
heat  flux  is  the  heat  flux  from  the  disturbance  not  including  the  additional  heat 
flux  from  the  electrical  heating  in  the  thin-film  sensor.  However,  it  is  readily 
apparent  that  at  best  the  apparatus  measures  an  average  heat  flux  into  the  front 
surface  of  the  gage  itself  in  the  region  of  the  thin  film.  Nevertheless,  if  the  average 
surface  temperature  measured  at  the  gage  approximates  the  surface  temperature 
that  would  exist  at  the  test  specimen  alone  and  the  electrical  heating  and  its 
effect  on  the  thermal  boundary  layer  is  small,  then  the  surface  heat  flux  values 
obtained  using  the  instrumented  test  specimen  approximate  the  external  surface 
heat  flux  that  would  exist  at  an  uninstrumented  test  specimen.  From  this  point 
on,  this  analysis  will  not  investigate  how  well  the  surface  heat  flux  measured  at  the 
gage  approximates  the  external  surface  heat  flux  that  would  be  present  at  the  test 
specimen  without  the  gage.  Instead,  the  analysis  investigates  how  well  the  series 
solution,  which  assumes  heat  transfer  in  an  initially  isothermal,  semi-infinite  solid, 
approximates  the  actual  external  surface  heat  flux  at  the  gage. 

1.2  Purpose 

The  purpose  of  this  investigation  is  to  numerically  simulate  the  heat  transfer 
in  a  thin-film  gage  with  a  finite-difference  model.  The  simulation  will  produce 
surface  temperatures  which  will  then  be  used  in  the  series  solution  to  find  estimates 
of  the  external  surface  heat  flux.  The  heat  fluxes  from  the  simulation  and  the  series 
solution  will  then  be  compared  to  evaluate  the  accuracy  of  the  series  solution. 

1.3  Approach 

Although  the  actual  heat  transfer  in  a  thin-film  gage  is  fully  three-dimen¬ 
sional,  the  analysis  uses  a  simplified  model  which  limits  the  heat  transfer  to  two 
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dimensions.  The  heat  transfer  is  limited  to  two  dimensions  by  modeling  the  region 
of  electrical  heat  generation  as  a  disk  centered  on  the  front  surface  of  the  quartz 
cylinder  rather  than  the  actual  thin  strip.  The  two  dimensions  in  the  heat  transfer 
problem  axe  now  the  radial  dimension  and  the  axial  dimension  into  the  quartz 
cylinder  from  the  front  surface.  Figure  2  illustrates  the  simplified,  two-dimensional 
geometry. 

The  series  solution  used  to  reduce  the  changes  in  surface  temperature  to 
values  for  surface  heat  flux  is  based  on  the  assumption  that  the  quartz  cylinder 
behaves  as  an  initially  isothermal,  semi-infinite  solid  when  the  temperature  data  is 
being  collected.  A  number  of  variables  influence  the  extent  to  which  the  cylinder 
departs  from  a  semi-infinite  solid.  These  variables  include  the  following: 

•  Geometry  parameters  such  as 

-  the  rati  of  the  length  of  the  quartz  cylinder  to  its  radius. 

-  the  area  covered  by  the  thin-film  sensor  compared  to  the  total  surface 
area  of  the  front  surface  of  the  quartz  cylinder.  (In  the  two-dimensional 
model,  the  ratio  of  the  radius  of  the  heated  disk  to  the  radius  of  the 
cylinder  describes  this  variable.) 

•  The  magnitude  of  the  thermal  disturbance. 

•  The  magnitude  of  the  electrical  heating  in  the  thin-film  sensor. 

•  The  actual  or  effective  thermal  properties  of  the  surrounding  materials. 

For  the  cylinder  to  behave  as  an  initially  isothermal,  semi-infinte  solid,  the 
temperature  distribution  in  the  cylinder  must  be  uniform  at  the  start  of  test  time, 
the  heat  transfer  in  the  cylinder  must  be  one-dimensional,  and  the  cylinder  must 
be  longer  than  the  distance  into  the  cylinder  to  which  the  effects  of  the  thermal 
disturbance  at  the  front  surface  propagate.  None  of  these  assumptions  may  be 
correct  in  an  actual  laboratory  experiment. 
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Figure  2. 


Simplified,  Two-Dimensional  Geometry 
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The  first  assumption,  that  the  initial  temperature  distribution  is  uniform,  is 
never  exactly  correct  because  the  electrical  heat  generation  in  the  thin-film  during 
the  time  when  the  bridge  is  being  balanced  establishes  a  non-uniform  temperature 
distribution.  The  magnitude  of  the  non-uniform,  initial  temperature  distribution 
depends  upon  the  magnitude  of  the  electrical  heat  generation,  the  external  condi¬ 
tions  and  the  amount  of  time  allowed  for  the  heat  generation  before  the  start  of 
test  time. 

The  second  assumption,  that  heat  transfer  is  one-dimensional,  is  also  never 
correct.  First  of  all,  the  initial  temperature  distribution  established  by  the  heat 
generation  before  the  start  of  the  test  is  multidimensional.  Therefore,  the  heat 
transfer  in  the  cylinder  departs  from  one- dimensional  from  the  very  start.  In 
the  simplified  model,  the  initial  temperature  distribution  will  be  two-dimensional. 
Secondly,  even  in  the  absence  of  preheating,  the  localized  surface  heat  generation 
and  the  fact  that  the  cylinder  is  not  ideally  isolated  from  the  surrounding  materials 
will  cause  the  one-dimensional  assumption  to  fail  eventually.  The  localized  surface 
heat  generation  will  produce  three-dimensional  heat  transfer  in  general,  but  the 
heat  transfer  is  limited  to  two  dimensions  by  the  simplified  geometry  of  the  model. 

The  third  assumption,  that  the  cylinder  is  long  compared  to  the  distance 
that  the  leading  edge  of  the  disturbance  travels,  is  only  valid  for  times  less  than 
the  time  it  takes  the  leading  edge  of  the  thermal  disturbance  to  travel  the  length 
of  the  quartz  cylinder.  A  rule  of  thumb  for  the  distance  traveled  by  the  leading 
edge  of  a  thermal  disturbance  into  a  solid  is  [l:pp. 60-61] 

/  »  4v/ot  (1) 

where  l  is  the  distance  to  which  the  leading  edge  of  the  disturbance  has  traveled, 
a  is  the  thermal  diffusivity  of  the  solid  and  t  is  the  time  after  exposure  to  the 
disturbance.  Using  this  rule  of  thumb,  one  would  expect  the  assumption  that  the 
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(2) 


cylinder  is  long  to  be  valid  only  for  test  times  less  than 

L2 

where  L  is  the  axial  length  of  the  cylinder. 

The  extent  to  which  the  cylinder  departs  from  a  semi-infite  solid  is  hard  to 
estimate.  It  is  even  more  difficult  to  estimate  the  error  in  the  values  for  exter¬ 
nal  surface  heat  flux  found  using  the  series  solution  with  the  surface  temperature 
changes  from  the  cylinder.  Because  it  is  difficult  to  analytically  estimate  the  mul¬ 
tidimensional  heat  transfer  in  the  cylinder  and  its  effect  on  the  results  for  external 
surface  heat  flux,  this  analysis  estimates  the  effects  through  simulation.  The  sim¬ 
plified,  two-dimensional  heat  transfer  problem  is  simulated  using  a  finite-difference 
model  on  a  digital  computer.  The  surface  temperatures  from  the  simulation  are 
used  to  find  the  series  solution  estimate  for  surface  heat  flux  which  is  then  com¬ 
pared  to  the  external  surface  heat  flux  of  the  simulation.  The  finite-difference 
model  allows  the  variables  to  be  changed  arbitrarily  so  that  one  can  investigate 
the  effect  of  each  variable  individually. 

In  investigating  the  effect  of  heat  transfer  across  the  outer-radial  boundary 
of  the  quartz  cylinder,  one  can  first  bracket  the  results  by  looking  at  three  limiting 
cases.  In  the  first  limiting  case,  one  can  model  the  cylinder  as  fully  insulated  at 
the  outer-radial  boundary  so  that  no  heat  flows  across  the  boundary.  This  limiting 
case  is  equivalent  to  surrounding  the  quartz  cylinder  by  a  material  whose  thermal 
conductivity  is  zero.  In  the  second  limiting  case,  one  can  surround  the  cylinder  with 
a  material  whose  thermal  conductivity  is  infinite.  The  second  material  is  assumed 
to  be  in  thermal  contact  with  a  heat  sink  maintained  at  the  initial  temperature 
of  the  gage.  In  this  case,  the  heat  flux  radially  out  of  the  cylinder  is  a  theoretical 
maximum.  In  the  third  limiting  case,  one  can  again  surround  the  cylinder  with 
a  material  whose  thermal  conductivity  is  infinite.  However,  for  this  case,  the 
temperature  of  the  surrounding  material  is  assumed  to  be  that  of  the  external 
disturbance.  In  this  case,  the  heat  flux  radially  into  the  cylinder  is  a  theoretical 
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maximum.  The  first  two  cases  bracket  all  possible  effects  due  to  heat  flux  radially 
out  of  the  cylinder.  The  first  and  third  cases  bracket  all  possible  effects  due  to 
heat  flux  radially  into  the  cylinder. 

In  addition  to  bracketing  the  results,  the  analysis  will  consider  a  couple  of 
intermediate  cases  for  heat  flux  across  the  outer-radial  boundary  of  the  cylinder. 
In  order  for  heat  to  flow  across  the  outer-radial  boundary  of  the  cylinder,  the  sur¬ 
rounding  material  must  have  a  non-zero  thermal  conductivity,  and  there  also  must 
be  a  temperature  gradient  at  the  boundary.  The  value  of  the  thermal  diffusivity 
in  the  surrounding  material  for  the  most  part  determines  the  direction  of  the  heat 
flux  established  across  the  outer-radial  boundary  of  the  cylinder.  Equation  1  shows 
that  a  thermal  disturbance  will  propagate  faster  in  a  material  with  a  larger  thermal 
diffusivity.  In  other  words,  a  material  with  a  larger  thermal  diffusivity  will  heat 
up  faster.  Thus,  if  the  thermal  diffusivity  of  the  surrounding  material  is  smaller 
than  the  thermal  diffusivity  of  the  quartz  cylinder,  heat  will  flow  radially  out  of  the 
cylinder.  On  the  other  hand,  if  the  thermal  diffusivity  of  the  surrounding  material 
is  greater  than  the  thermal  diffusivity  of  the  quartz  cylinder,  heat  may  flow  into 
the  cylinder.  It  is  important  to  note,  however,  that  the  electrical  heat  generation 
will  always  tend  to  establish  a  temperature  gradient  for  heat  to  flow  radially  out¬ 
ward.  For  heat  to  flow  into  the  cylinder,  this  temperature  gradient  must  first  be 
overcome. 

Another  simplification  in  the  two-dimensional  model  involves  the  way  in 
which  the  intermediate  cases  for  heat  flux  across  the  outer-radial  boundary  of 
the  cylinder  are  modeled.  As  mentioned  in  the  first  section,  the  gage  is  embedded 
in  the  test  specimen  using  a  filler  material.  Thus,  the  true  heat  transfer  problem 
involves  three  materials:  the  quartz  cylinder,  the  filler  material,  and  the  test  spec¬ 
imen.  However,  the  surrounding  materials  affect  the  results  for  external  surface 
heat  flux  only  to  the  extent  that  they  induce  heat  flux  into  or  out  of  the  cylinder. 
The  model  for  the  intermediate  cases  for  heat  flux  across  the  outer-radial  bound- 
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ary  of  the  cylinder  uses  only  two  materials:  the  quartz  cylinder  and  a  surrounding 
material.  Using  two  materials  is  equivalent  to  modeling  the  effective  tendency  of 
the  surrounding  materials  to  induce  heat  flux  across  the  outer-radial  boundary  of 
the  cylinder. 
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II.  The  Finite- Difference  Model 


Before  discussing  the  finite-difference  model,  it  is  instructive  to  first  look  at 
the  differential  equations  for  the  two-dimensional  model.  A  good  understanding 
of  the  differential  equations  will  help  in  developing  the  finite-difference  model  and 
in  understanding  the  approximations  caused  by  the  discretization.  Also,  non- 
dimensionalizing  the  differential  equations  will  yield  the  important  non-dimensional 
parameters  in  the  two-dimensional  model. 

2.1  The  Differential  Equations 

Incropera  and  Dewitt  derive  the  governing  equations  for  heat  transfer  in 
a  solid  using  an  energy  balance  on  the  appropriate  differential  control  volume 
[4:pp. 43-53].  For  a  homogeneous,  isotropic  solid — one  in  which  properties  are 
constant  with  position  and  uniform  in  all  directions— the  governing  equation  for 
temperature  in  two-dimensional,  cylindrical  coordinates  is 

id_(  ot\  |  o2t  _  _iar 

r  dr  y  dr)  dz 2  a  dt 

where  r  is  the  radial  dimension,  z  is  the  axial  dimension  and  thermal  diffusivity, 
a,  is  the  ratio  of  thermal  conductivity  to  thermal  capacitance,  a  =  k/pcp  .  A 
general  solution  of  Equation  3  will  give  the  temperature  distribution  in  the  solid 
as  a  function  of  the  two  spacial  coordinates,  z  and  r,  and  time  t.  The  solution  of 
Equation  3  is  governed  by  the  initial  condition  on  time  and  the  boundary  conditions 
for  each  of  the  two  spacial  coordinates. 

The  initial  condition  for  the  two-dimensional  model  before  preheating  is  sim¬ 
ply 

T(r,  z)  =  T,  for  t  <  0  (4) 

The  boundary  conditions  for  the  two-dimensional  model  require  further  explana¬ 
tion. 
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2.1.1  Front  Surface  Boundary  Condition  The  front  surface  boundary  con¬ 
dition  must  include  consideration  of  the  electrical  heat  generation  in  the  thin-film 
sensor.  As  previously  mentioned,  the  two-dimensional  model  models  the  sensor 
as  a  disk  centered  on  the  front  surface  of  the  cylinder.  Using  three  assumptions, 
the  thin-film  sense  can  be  modeled  simply  as  a  region  of  surface  heat  generation. 
First,  because  the  sensor  is  an  extremely  thin  film  of  platinum,  the  total  thermal 
capacity  of  the  sensor  is  small.  Secondly,  the  thermal  contact  between  the  plat¬ 
inum  film  and  the  quartz  cylinder  is  assumed  to  be  good  so  that  the  temperature 
of  the  lower  side  of  the  platinum  film  is  the  same  as  the  surface  temperature  of  the 
quartz  below  it.  Lastly,  the  thermal  conductivity  of  platinum  is  large.  Because  the 
sensor  is  very  thin  and  has  a  large  thermal  conductivity,  the  temperature  through 
the  sensor  essentially  will  be  constant  and  equal  to  the  temperature  of  the  quartz 
below  it.  Because  the  thermal  capacity  of  the  sensor  is  small,  its  energy  storage  is 
negligible.  Figure  3  illustrates  an  energy  balance  across  the  sensor  using  the  above 
assumptions.  Performing  the  energy  balance  yields 
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where  q,  is  the  heat  flux  into  the  cylinder,  hFRONT{TjFROf]J.  —  T, )  is  the  heat  flux 
into  the  disk  from  the  fluid  and  qg  is  the  electrical  heat  generation  per  surface  area 
of  the  disk.  Then,  the  boundary  condition  for  the  front  surface  of  the  cylinder  is 
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2.1.2  Back  Surface  Boundary  Condition  The  boundary  condition  at  the 
back  surface  is  similar  to  the  boundary  condition  at  the  front  surface  without  the 
complication  of  the  thin-film  sensor.  Us:ag  the  sign  convention  that  positive  heat 
flux  points  in  the  positive  z  direction,  positive  heat  flux  at  the  back  surface  is  heat 
flux  out  of  the  cylinder  instead  of  into  the  cylinder: 

96  =  ~k  ~ KacMback  ~  Tfc)  (7) 
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Figure  3.  Energy  Balance  Across  the  Thin-Film  Sensor 


where  T, j  is  the  temperature  on  the  back  surface  of  the  cylinder. 


£.1.8  Inner- Radial  Boundary  Condition  For  the  two-dimensional  model, 
heat  flux  in  the  radial  direction  anywhere  in  the  cylinder  is  the  same  regardless  of 
the  angular  position  in  the  cylinder.  This  radial  symmetry  can  be  satisfied  only  if 
the  heat  flux  in  the  radial  direction  at  the  centerline  is  zero.  Then,  the  inner- radial 
boundary  condition  for  the  two-dimensional  model  is 


9r 


(8) 


2.1.4  Outer- Radial  Boundary  Condition  A  general  statement  for  the  bound¬ 
ary  condition  at  the  outer  radius  of  the  cylinder  is  similar  to  the  general  statement 
for  the  boundary  condition  at  the  front  and  back  surfaces.  The  heat  flux  crossing 
the  outer-radial  boundary  must  be  equal  to  the  heat  flux  in  the  radial  direction  on 
either  side  of  the  boundary.  In  this  case,  however,  both  sides  of  the  boundary  are 
solids,  so  the  boundary  condition  uses  Fourier’s  law  for  heat  conduction  on  both 
sides: 


9r 


-t  -  -*■ 


r-RcyL 


(9) 


dr )  R_  dr 

where  k'  is  the  thermal  conductivity  of  the  surrounding  material.  An  important 
assumption  in  using  this  boundary  condition  is  that  the  thermal  contact  between 
the  quartz  cylinder  and  the  surrounding  material  is  good. 


As  described  in  Section  1.3,  it  is  worthwhile  to  investigate  three  limiting  cases 
for  the  boundary  condition  at  the  outer  radius  of  the  cylinder.  Each  of  the  limiting 
cases  can  be  modeled  without  actually  including  the  surrounding  material  in  the 
model. 


The  first  limiting  case  is  the  case  for  no  heat  flux  across  the  outer-radial 
boundary  of  the  cylinder.  The  boundary  condition  for  the  fully-insulated  or  adia¬ 
batic  case  is 


?r 


—  k 


ar\ 

dr)  R 

/  r=RCYL 


=  0 


(10) 
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The  second  limiting  case  is  the  case  for  theoretically  maximum  heat  flux  out 
across  the  outer-radial  boundary  of  the  cylinder.  This  limiting  case  is  equivalent  to 
assuming  that  the  thermal  conductivity  of  the  surrounding  material  is  infinite  and 
that  the  surrounding  material  is  in  thermal  contact  with  a  heat  sink  maintained 
at  the  initial  temperature  of  the  gage.  Then,  the  temperature  in  the  surrounding 
material  is  always  the  initial  temperature.  Thus,  the  outer-radial  boundary  con¬ 
dition  on  the  cylinder  for  the  limiting  case  of  maximum  heat  flux  out  is  simply  a 
condition  of  constant  temperature: 

T(r  =  RCYL,z)  =  T,  (11) 

The  third  limiting  case  is  the  limiting  case  for  theoretically  maximum  heat 
flux  in  across  the  outer-radial  boundary.  This  case  is  equivalent  to  assuming  that 
the  thermal  conductivity  of  the  surrounding  material  is  infinite  and  that  the  sur¬ 
rounding  material  is  maintained  at  the  temperature  of  the  external  disturbance. 
The  outer-radial  boundary  condition  on  the  cylinder  for  the  limiting  case  of  max¬ 
imum  heat  flux  in  is  also  a  condition  of  constant  temperature: 

T(r  =  RCYL,z)  =  T}  (12) 

When  including  the  surrounding  material  in  the  model,  the  differential  equa¬ 
tions  for  heat  transfer  in  the  surrounding  material  are  identical  to  those  for  heat 
transfer  in  the  cylinder  except  that  the  properties  for  the  surrounding  material  are 
used  in  place  of  the  properties  for  quartz.  The  outer-radial  boundary  for  the  sys¬ 
tem  is  now  the  outer-radial  boundary  of  the  surrounding  material,  and  one  of  the 
three  limiting  cases  must  be  used  for  the  outer-radial  boundary  condition  on  the 
system.  Equation  9  then  serves  as  a  compatability  equation  at  the  discontinuity 
between  the  two  materials. 
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2.2  Non-Dimensional  Parameters 


When  non-dimensionalizing  the  differential  equations,  it  is  first  necessary  to 
specify  the  non-dimensional  form  for  the  variables.  Then,  the  form  of  the  equa¬ 
tions  will  reveal  the  relevant  non-dimensional  parameters  of  the  problem.  In  this 
problem,  it  is  helpful  to  separate  the  total  heat  transfer  problem  into  two  sepa¬ 
rate  problems  which  use  different  non-dimensional  forms  for  the  variables  and  have 
some  different  non-dimensional  parameters. 

The  first  of  the  two  problems  is  the  preheating  problem.  The  preheating 
problem  models  the  period  of  time  during  which  the  instrumentation  circuitry  is 
turned  on  and  the  bridge  is  being  balanced.  The  test  specimen  is  not  exposed  to 
the  external  disturbance  during  the  preheating  problem.  Instead,  the  system  is 
disturbed  by  the  electrical  heat  generation.  There  are  three  basic  assumptions  for 
the  preheating  problem: 

•  The  external  fluid  temperature  on  both  the  front  and  back  surfaces  remain 
at  the  initial  temperature  of  the  system. 

•  The  values  of  the  convection  coefficients  on  both  surfaces  are  typical  free 
convection  values. 

•  The  surface  heat  generation  is  always  some  non-zero  value. 


For  the  preheating  problem,  the  non-dimensional  forms  for  the  two  dependent 
variables  T  and  q  are 


T-Tj 

i<lg/p  h FRONT ) 

s_ 

<19 


(13) 

(14) 


The  second  problem  is  the  disturbance  problem.  The  disturbance  problem 
models  the  heat  transfer  after  the  test  specimen  is  exposed  to  the  external  dis¬ 
turbance.  The  initial  temperature  distribution  for  the  disturbance  problem  is  the 
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final  temperature  distribution  from  the  preheating  problem.  There  are  three  basic 
assumptions  for  the  disturbance  problem  as  well: 


•  The  external  fluid  temperature  at  the  front  of  the  test  specimen  is  different 
from  the  initial  temperature  of  the  system. 

•  The  value  of  the  convection  coefficient  on  the  front  surface  of  the  test  speci¬ 
men  can  be  any  typical  forced  convection  value. 

•  Heat  generation  may  or  may  not  be  present. 


For  the  disturbance  problem,  the  non-dimensional  forms  for  the  two  dependent 
variables  are 


d0+ 


dq 


+ 


T-Ti 


dTj  FRONT  T 

q 


d'1  FRONT 


UT,t 


FRONT 


-Ti) 


(15) 

(16) 


It  will  be  valuable  to  look  at  cases  where  the  heat  generation  is  zero.  If  the 
surface  heat  generation  is  zero,  then  the  solution  to  the  preheating  problem  is  the 
trivial  solution;  the  temperature  remains  at  the  initial  temperature  of  the  system 
until  exposed  to  the  external  disturbance.  Thus,  the  total  heat  transfer  problem 
includes  only  the  disturbance  problem  when  the  heat  generation  is  zero. 

The  preheating  problem  and  the  disturbance  problem  are  non-dimensional- 
ized  differently  to  account  for  the  different  conditions  in  the  two  problems.  The 
preheating  problem  is  non-dimensionalized  with  qg/ dhFKOST  denomina¬ 

tor  because  the  temperature  difference  dTfFRONT  —  T  is  zero  in  the  preheaing 
problem  but  the  heat  generation  is  never  zero.  The  disturbance  problem  is  non- 
dimensionalized  with  the  temperature  difference  dT/FRONT  -  T,  in  the  denominator 
instead  to  allow  the  surface  heat  generation  to  be  zero  for  the  disturbance  problem. 
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The  non-dimensional  forms  for  the  independent  variables  r,  z  and  t  are  the 


same  in  both  problems.  They  are 

r+  =  / 

^ CYL 

(17) 

*♦  -  i 

(18) 

-  % 

(19) 

With  these  definitions  for  the  non-dimensional  variables,  the  differential  equa¬ 
tions  can  be  non-dimensionalized.  The  non-dimensional,  governing  equation  for 
temperature  within  the  solid  for  either  problem  is 


[J 

lA 


1  2 


d29+  d9+ 

r+  dr+  \^r  dr+  )  dz+ 2  dt+ 


(20) 


The  initial  condition  for  the  preheating  problem  is 


0+(r+,  2+)  =  0  for  t  <  0 


(21) 


Recall  that  the  initial  condition  for  the  disturbance  problem  is  the  final  temperature 
distribution  from  the  preheating  problem. 

It  is  interesting  to  note  that  Equations  (20)  and  (21)  remain  the  same  regard¬ 
less  of  the  non-dimensional  form  for  temperature  difference,  9  =  T  —  T)  .  One  can 
multiply  Equations  (20)  and  (21)  by  my  non-dimensional  constant  changing  the 
form  for  the  non-dimensional  variable,  9 ,  and  the  equations  will  remain  the  same. 

The  non-dimensional  equations  for  the  two  seperate  problems  differ  only  in 
some  of  the  boundary  conditions.  The  non-dimensional  boundary  conditions  for 
the  preheating  problem  are  as  follows: 


Top  surface  boundary  condition; 


P^FRONT^ 


1  -1 


d9+' 

dz+ 


z+=0 


|  -6+  +  1  for  0  <  r+  < 

-9t 


for  [tPH  <  r+  < 

KcyL  —  KcYL 


(22) 
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Bottom  surface  boundary  condition; 


phBACKL]~l  dO+ 


dz+J^=x 


Inner-radial  boundary  condition; 


d$+\ 

*+A+-o~ 


Outer-radial  boundary  conditions ; 


Boundary  between  materials  — for  use  when  including  the  surrounding  ma¬ 
terial  in  the  model, 

de+\  rjt'l  d0+\  ,  s 

dr+)T+=i-  ~  M  dr+)r+=i+  25 


Fully-insulated, 


Maximum  heat  flux  out, 


Maximum  heat  flux  in, 


for  the  limiting  cases. 


dr+)  r+=l 


e+(r+  =  l,z+)  =  0 


e+(r+  =  l,z+)  =  0 


The  non-dimensional  boundary  conditions  for  the  disturbance  problem  are 
the  following: 


Top  surface  boundary  condition; 


dhFRONTL\  d6+\ 

k  dz+)  ^=Q 


1  -  Of  + 

i-e 


(' laldhfKQHT ) 

Off  pRont~T') 


for  0  <  r+  <  ] 

fOT  tel  <  ^  s  tef] 


Bottom  surface  boundary  condition; 


{  b’  n  .  r*  *.  L 


1  -1 


ajL\  -«♦ 


(31) 


Inner-radial  boundary  condition; 


dO+ 

dr+ 


^  r+=0 


(32) 


Outer-radial  boundary  conditions; 


Boundary  between  materials  — for  use  when  including  the  surrounding  ma¬ 
terial  in  the  model, 


do+\ 

dr+Jr+  =  l- 


d#+\ 

dr+/r+=  1  + 


(33) 


or 


Fully-insulated, 


Maximum  heat  flux  out, 


Maximum  heat  flux  in 


for  the  limiting  cases. 


d0+\ 


=  0 


5r+/r+  =  , 

e+(r+  =  1,2+)  =  0 

0+(r+  =  1  ,z+)  =  1 


(34) 


(35) 


(36) 


The  non-dimensional  differential  equations  for  the  heat  transfer  in  the  cylin¬ 
der  during  the  preheating  and  disturbance  problems  include  the  following  non- 
dimensional  parameters  which  govern  the  solution: 

External  fluid  conditions  and  electrical  heat  generation; 


pBifRONT 


phFRONTL 
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(37) 
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P^lBACK 

P^BACK^ 

k 

(38) 

front 

d^FRONT^ 

k 

(39) 

BACK 

d^BACK^ 

k 

(40) 

7 

(?S / dh  front) 

(dTfFRON T  ~  Ti) 

(41) 

T 

( dTfBACK  -  T{) 

( dT!  FRONT  ~  Ti) 

(42) 

Geometry  ratios; 

L 

T> 

CYL 

(43) 

n 

IXD1SK 

D 

Acyi 

(44) 

Property  ratio 

k' 

k 

(45) 

Two  additional  non 

-dimensional  parameters  arise 

when  including  the  sur- 

rounding  material  in  the  analysis.  The  first  is  a  geometry  ratio  specifying  the 
extent  of  the  surrounding  material: 


kmax 


R 


CYL 


The  second  is  another  property  ratio — the  ratio  of  thermal  capacitance 


(pep) 

(pcvY 


(46) 


(47) 


When  incorporating  the  final  temperature  distribution  from  the  preheating 
problem  as  the  initial  temperature  distribution  to  the  disturbance  problem,  an¬ 
other  non-dimensional  parameter  can  be  defined  for  convenience.  Non-dimensional 
temperature  values  from  the  preheating  problem  can  be  converted  to  the  non- 
dimensional  form  used  in  the  disturbance  problem  by  multiplying  by  the  generation 
ratio,  7,  and  a  derived  parameter,  /?,  which  is  the  ratio  of  the  convection  coefficient 
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at  the  front  surface  in  the  disturbance  problem  to  the  convection  coefficient  at  the 
front  surface  in  the  preheating  problem: 

P  =  (48) 

FRONT 

2.3  The  Finite- Difference  Equations 

Without  an  analytical  solution  to  the  differential  equations  for  the  two- 
dimensional  model,  it  is  necessary  to  approximate  the  transient  solution  through 
simulation.  The  most  common  method  for  simulating  heat  transfer  is  to  run  a 
finite-difference  model  on  a  digital  computer.  Various  finite-different  schemes  are 
available.  This  analysis  uses  an  explicit  finite-difference  scheme  derived  using  the 
control- volume  approach.  The  explicit  scheme  uses  second  order,  central  differences 
for  the  spatial  derivatives  and  first  order,  forward  differences  for  the  time  deriva¬ 
tives.  Incropera  and  Dewitt  present  the  control-volume  approach  in  developing 
finite-difference  models  for  conduction  heat  transfer  problems  [4:pp.  143-149, 213- 
220].  The  control-volume  approach  has  the  advantages  of  being  easy  to  use  and 
versatile  in  application  to  different  and  complicated  boundary  and  geometry  con¬ 
ditions.  The  control-volume  approach  is  also  a  conservative  approach  which  means 
that  the  finite-difference  model  developed  using  the  control-volume  approach  will 
satisfy  the  statement  of  conservation  of  energy  to  within  the  truncation  error  for 
the  model.  The  finite-difference  models  used  for  the  single-material  and  the  two- 
material  problems  with  fully-insulated,  outer-radial  boundary  conditions  are  given 
in  Appendices  D.l  and  D.2,  respectively. 

2.3.1  Discretizing  the  Model  The  first  step  in  developing  the  finite-differ¬ 
ence  model  is  to  discretize  the  continuous  model  in  both  space  and  time.  Figures  4 
and  5  illustrate  the  discretized  geometry  of  the  single-material  problem  and  the 
two-material  problem,  respectively.  To  discretize  the  geometry,  the  cylinder  and 
the  surrounding  material  are  subdivided  into  nodal  regions.  Associated  with  each 
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Figure  5.  Discretized  Geometry  of  the  Two-Material  Problem 
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nodal  region  is  a  reference  point  or  reference  level  for  the  region  called  a  node.  All 
nodal  regions  are  annular  rings  except  for  those  on  the  centerline  of  the  cylinder 
which  are  cylinders  of  radius  Ar/2.  Interior  nodal  regions — those  not  on  an  exterior 
system  boundary — are  annualar  rings  of  width  Ar  and  depth  A 2.  Nodal  regions  on 
exterior  system  boundaries  are  generally  assigned  half  the  step  size  of  interior  nodal 
regions.  For  nodal  regions  at  any  boundary,  the  node  is  placed  on  the  boundary. 
Placing  nodes  directly  on  the  boundary  rather  than  one-half  spacial  step  in  from 
the  boundary  provides  better  resolution  of  the  boundary  conditions.  Because  the 
heat  transfer  in  the  two-dimensional  model  is  axi-symmetric,  the  reference  point 
for  nodal  regions  other  than  those  on  the  centerline  of  the  cylinder  can  be  any 
point  on  an  annular  ring  at  the  reference  radius  and  depth. 

In  general,  discretizing  the  geometry  means  that  the  finite-difference  model 
only  solves  for  temperatures  at  the  nodes.  The  temperature  at  a  node  represents 
an  average  temperature  for  the  nodal  region.  For  the  purpose  of  deriving  the 
finite-difference  equations,  conditions  across  the  nodal  region  are  often  assumed 
constant  and  equal  to  the  conditions  associated  with  the  node  or  the  nodal  region. 
Continuous  time  is  discretized  similarly,  so  the  finite-difference  model  only  solves 
for  temperatures  at  discrete  intervals  in  time  as  well. 

To  identify  positions  in  the  nodal  geometry  mesh,  the  index  value  n  specifies 
the  radial  position  of  the  node,  and  the  index  value  m  specifies  the  axial  position  of 
the  node.  When  used  as  subscripts,  the  radial  index  is  placed  first  such  as  in  0n  m. 
(The  only  exceptions  are  for  values  denoting  surface  conditions  using  the  subscripts 
*  and  6  to  denote  the  axial  positions  on  the  front  and  back  surfaces,  respectively. 
Subscripts  3  and  b  are  placed  before  the  radial  index  value  as  in  #,,„•)  The  index 
value  of  zero  for  n  denotes  postions  on  the  centerline  of  the  cylinder  where  the 
radial  dimension  is  zero.  The  index  value  of  zero  for  m  denotes  positions  on  the 
front  surface  where  the  axial  dimension  is  zero.  Then,  the  radial  position  of  a  node 
with  radial  index  n  is  r  =  nAr  .  Similarly,  the  axial  position  of  a  node  with  axial 
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index  m  is  z  —  mAz  .  The  index  value  j  specifies  the  time  level  in  discretized 
time.  The  time  index  j  is  usually  placed  as  a  superscript.  The  index  value  of  zero 
for  j  denotes  time  t  =  0  ,  so  the  *ime  denoted  by  index  j  is  t  =  j  At  . 

2.8.2  Nodal  Equations  In  the  explicit  scheme,  a  finite-difference  equation 
for  each  node  estimates  the  temperature  of  the  node  at  the  new,  unknown  time 
level  in  terms  of  the  temperatures  of  the  node  itself  and  surrounding  nodes  at 
the  old,  known  time  level.  The  finite-difference  equations  are  developed  using  an 
energy  balance  on  the  individual  nodal  volumes. 

For  the  purpose  of  deriving  the  finite- difference  equations,  one  assumes  that 
the  heat  flux  at  each  surface  of  the  nodal  volume  is  directed  into  the  nodal  volume. 
Of  course,  the  actual  direction  of  the  heat  flux  during  the  simulation  depends 
on  the  temperature  profile.  This  sign  convention  merely  helps  in  developing  the 
equations  correctly.  One  also  assumes  that  the  heat  flux  at  each  surface  of  the 
nodal  volume  is  constant  and  equal  to  the  heat  flux  at  the  reference  level.  Then, 
the  energy  balance  at  each  nodal  volume  is 

q\Ai  +  ^2^2  +  <73^3  +  <?4-44  =  (Pcp)V  (49) 

The  subscripts  1-4  denote  the  inner-radial,  outer-radial,  front  (smaller  axial  dimen¬ 
sion)  and  back  (larger  axial  dimenion)  surfaces  of  the  nodal  volume,  respectively. 

The  develoment  of  each  nodal,  finite-difference  equation  requires  specific  in¬ 
formation  about  the  node  whose  equation  is  being  developed.  Nevertheless,  the 
procedure  is  basically  the  same  for  all  nodes.  To  illustrate  the  process,  the  steps  are 
outlined  for  an  interior  node  in  the  cylinder.  Figure  6  shows  the  nodal  geometry 
for  an  interior  node  in  the  cylinder. 

The  first  step  is  to  specify  the  dimensions  of  the  nodal  region.  For  an  interior 
node  in  the  cylinder 

A\  =  2tt  (n  —  l/2)ArAz  (50) 
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Ai  =  2?r(n  +  l/2)ArAz  (51) 

Az  =  7r  |(nAr  +  Ar/2)2  -  (nAr  -  Ar/2)2]  =  7r(2n)Ar2  (52) 

A4  =  7r(2n)Ar2  (53) 

and 

V  =  7r(2n)Ar2Az  (54) 


The  next  step  is  to  determine  a  finite-difference  expression  for  the  heat  flux 
at  each  of  the  surfaces  of  the  nodal  volume.  For  an  interior  node,  a  finite-difference 
expression  for  the  heat  flux  at  each  surface  is  determined  from  Fourier’s  conduction 
law  using  a  central  difference  to  approximate  the  temperature  gradient.  Using  the 
finite-difference  expressions  for  heat  flux  and  a  forward  difference  to  approximate 
the  time  derivative,  a  finite-difference  expression  for  the  energy  balance  at  an 
interior  node  in  the  cylinder  is 

r  n  a  2  a  1  (Tti-TiJ 
pCp  [7r(2n)Ar'Azj  - At~  '  = 

(Tj  ,  -T3  ) 

k  [2n (n  -  l/2)ArAz]  ’^r  — 

+  k  [27r(n  +  l/2)ArAz]  ^Tn+1'^~  (55) 

-|-  k  |7r(2n)Ar  j  - 1 1 —  +  k  [7r(2n)Ar  j  - — - 

By  calculating  the  spacial  derivatives  at  the  old  time  level,  j,  and  not  at 
the  new  time  level,  j  +  1,  or  some  combination  of  the  two,  the  finite-difference 
scheme  will  be  an  explicit  scheme.  The  only  unknown  in  each  nodal  equation  is 
the  temperature  of  the  node  at  the  new  time  level,  j  +  1,  used  in  the  forward- 
difference  approximation  for  the  time  derivative.  Solving  for  the  temperature  of 
each  node  at  the  new  time  level  is  straightforward.  Grouping  temperature  terms 
for  the  same  node  at  the  same  time  level  toget'  er  and  solving  for  the  temperature 
at  the  new  time  level  gives  the  dimensional  form  for  the  nodal  finite-difference 


equation  for  an  interior  node  in  the  cylinder: 


0(1  -  +°<i + 

At  At 

+  a(Al)2Tn,m“1  +  a(A^T"’m+1 


(56) 


It  is  not  necessary  to  derive  the  finite-difference  equation  at  each  node.  Nodes 
that  have  similar  geometry  and  the  same  boundary  conditions  will  have  identical 
finite-difference  equations.  The  only  difference  in  the  equations  will  be  the  specific 
values  of  the  indexing  variables.  Also,  when  including  the  surrounding  material  in 
the  model,  one  must  be  careful  to  use  the  correct  property  values  for  each  term 
in  the  energy  balance.  For  the  single-material  problem,  there  are  nine  different 
forms  for  the  nodal  finite-difference  equations  corresponding  to  the  nine  different 
locations  labeled  (a)-(i)  in  Figure  4.  For  the  two-material  problem,  in  which 
both  the  cylinder  and  the  surrounding  material  are  included  in  the  finite-difference 
nodal  mesh,  there  are  fifteen  different  forms  for  the  nodal  finite-difference  equations 
corresponding  to  the  locations  labeled  (a)-(o)  in  Figure  5. 


2.3.3  Boundary  Conditions  An  advantage  of  the  control-volume  approach 
in  deriving  the  finite-difference  equations  is  the  ease  with  which  the  boundary  con¬ 
ditions  and  material  boundaries  are  incorporated  in  the  finite-difference  model. 
The  continuous  forms  of  the  boundary  conditions  for  the  two-dimensional  model 
were  introduced  in  Sections  2. 1.1-2. 1.4.  Recall  that  the  single-material  problem 
uses  one  of  the  three  limiting  cases  for  outer-radial  boundary  condition  on  the  cylin¬ 
der.  The  two-material  problem  must  also  use  one  of  the  three  limiting  cases  for 
the  outer-radial  boundary  condition  on  the  surrounding  material.  The  boundary 
between  the  materials  in  the  two-material  problem  will  require  special  considera¬ 
tion. 

To  incorporate  the  front  surface  boundary  condition  in  the  finite-difference 
model,  the  nodal  energy  balance  for  nodes  at  the  front  surface  is  modified.  The 
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heat  flux  into  the  nodal  region  across  the  front  surface,  q3,  is  determined  from  the 
front  surface  boundary  condition,  Equation  (6): 
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(57) 


The  only  approximation  in  applying  the  front  surface  boundary  condition  to  the 
finite-difference  model  is  the  assumption  that  the  heat  flux  is  constant  across  the 
nodal  region  and  equal  to  the  heat  flux  evaluated  at  the  node.  Assuming  that  the 
heat  flux  is  constant  is  equivalent  to  assuming  that  the  temperature  on  the  surface 
of  the  nodal  region  is  constant  and  equal  to  the  temperature  at  the  node. 


The  back  surface  boundary  conditon  is  incorporated  similarly.  The  nodal 
energy  balance  for  nodes  on  the  back  surface  uses  the  back  surface  boundary  con¬ 
dition  of  Equation  (7)  to  determine  the  heat  flux  into  the  node  across  the  back 
surface: 

q4  —  h  BACK{T jBACK  —  mmax)  (58) 


One  of  the  three  limiting  cases  for  the  outer-radial  boundary  condition  must 
be  used  for  the  outer-radial  boundary  condition  of  the  system.  To  use  the  fully- 
insulated  case,  the  energy  balance  for  nodes  on  the  outer-radial  boundary  is  mod¬ 
ified  by  setting  the  heat  flux  into  the  node  from  a  larger  radial  position  equal  to 
zero: 

<72  =  0  (59) 

The  cases  for  maximum  heat  flux  out  across  the  boundary  and  maximum  heat  flux 
in  across  the  boundary  are  realized  by  keeping  the  temperature  of  the  node  on  the 
boundary  equal  to  some  constant.  The  nodes  on  the  outer-radial  boundary  do  not 
need  a  finite-difference  equation  since  the  temperature  at  the  next  time  level  does 
not  change.  In  the  case  for  maximum  heat  flux  out,  the  temperature  at  each  node 
on  the  outer-radial  boundary  remains  at  the  initial  temperature: 


=  Ti 


(60) 
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In  the  case  for  maximum  heat  flux  in  across  the  outer-radial  boundary,  the  temper¬ 
ature  at  each  node  on  the  outer-radial  boundary  remains  at  the  fluid  temperature: 

=  2/  (61) 


2.3.4  The  Boundary  Between  Materials  For  the  two-material  problem,  in 
which  the  thermal  conductivities  of  the  materials  are  different,  the  temperature 
gradient  at  the  boundary  between  the  materials  is  discontinuous.  Nevertheless, 
energy  is  still  conserved.  Because  the  control- volume  approach  is  based  on  the  con¬ 
servation  of  energy  for  each  nodal  volume,  the  development  of  the  finite-difference 
equations  for  nodes  on  the  boundary  between  the  materials  remains  essentially  the 
same.  The  same  assumptions  are  used  with  special  care  to  account  for  the  change 
in  properties.  Figure  7  shows  the  nodal  geometry  for  a  node  on  the  boundary 
between  the  materials. 


The  conductivities  to  use  in  the  expression  for  heat  fluxes  qi  and  q2  are  readily 
apparent.  Heat  flux  qx  in  the  cylinder  uses  the  thermal  conductivity  of  the  cylinder, 
while  heat  flux  q2  in  the  surrounding  material  uses  the  thermal  conductivity  of  the 


surrounding  material: 
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The  heat  flux  at  the  front  and  back  surfaces,  <73  and  <?4,  is  a  function  of  both 
thermal  conductivities.  The  axial  temperature  gradient  at  each  surface  is  still 
assumed  constant  and  equal  to  the  axial  temperature  gradient  evaluated  at  the 
boundary.  However,  the  heat  flux  at  the  portion  of  the  surface  in  the  cylinder  is 
determined  using  the  thermal  conductivity  of  the  cylinder,  while  the  heat  flux  at  the 
portion  of  the  surface  in  the  surrounding  material  is  determined  using  the  thermal 
conductivity  of  the  surrounding  material.  The  heat  flow  through  each  portion  of 
each  surface  is  the  product  of  the  heat  flux  at  each  portion  of  each  surface  times 
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Figure  7.  Nodal  Geometry  for  a  Node  on  the  Boundary  Between  Materials 


the  area  of  the  portion  of  the  surface.  Then,  the  total  heat  flow  across  each  surface 
is  equal  to  the  sum  of  the  heat  flows  across  each  portion.  Finally,  the  average  heat 
flux  across  each  surface  is  the  total  heat  flow  divided  by  the  total  surface  area. 


An  equivalent  method  of  calculating  the  heat  flux  at  each  surface  is  to  use 
an  average  value  for  thermal  conductivity,  k,  which  is  weighted  according  to  the 
fraction  of  the  surface  area  in  each  of  the  two  materials: 


£  _  k(NCYL  —  1/4)  +  k'(NcYL  +  1/4) 
2  Ncyl 

Then,  the  finite-difference  expressions  for  93  and  94  are 
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The  total  thermal  capacity  of  the  nodal  volume  is  equal  to  the  thermal  capac¬ 
ity  contributed  by  the  surrounding  material  plus  the  thermal  capacity  contributed 
by  the  quartz  material.  The  thermal  capacity  contributed  by  the  surrounding 
material  is  the  thermal  capacitance  of  the  surrounding  material  times  the  part  of 
the  nodal  volume  made  up  by  the  surrounding  material.  The  thermal  capacity 
contributed  by  the  quartz  material  is  determined  similarly.  Again,  an  equivalent 
method  of  calculating  the  total  thermal  capacity  is  to  multiply  the  total  nodal  vol¬ 
ume  by  an  average  thermal  capacitance  which  this  time  is  weighted  by  the  fraction 
of  volume  in  each  material: 


Op)  = 
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(67) 


Then,  the  energ y  storage  term  for  nodal  volumes  on  the  boundary  between  the 
materials  is 
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2.8.5  Non-Dimensional  Nodal  Equations  The  nodal  equations  are  non-di- 
mensionalized  the  same  way  that  the  differential  equations  are  non-dimensionalized. 
First,  the  non-dimensional  form  of  the  variables  is  specified.  Then,  the  non- 
dimensional  form  of  the  equations  determines  the  relevant  non-dimensional  param¬ 
eters  in  the  problem.  As  explained  in  Section  2.2,  the  total  heat  transfer  problem 
is  divided  into  two  seperate  problems,  the  preheating  problem  and  the  disturbance 
problem.  The  two  problems  use  different  non-dimensional  forms  for  the  dependent 
variables  T  and  q  but  the  same  non-dimensional  forms  for  the  independent  vari¬ 
ables  r,  z  and  t.  The  non-dimensional  nodal  equations  are  the  same  except  for  the 
nodes  on  both  the  front  and  back  surfaces  (see  Section  2.2). 

The  non-dimensional  forms  for  temperature  and  heat  flux  in  the  preheating 
and  disturbance  problems  are  given  in  Equations  (13)-(16)  in  Section  2.2.  As 
in  the  dimensional  form  of  the  nodal  equations,  the  independent  variables  in  the 
non-dimensional  nodal  equations  are  expressed  as  an  index  value  multiplied  by  a 
non-dimensional  step  size: 


r+  =  nAr+  (69) 

2+  =  mAz+  (70) 

t+  =  jAt+  (71) 


The  non-dimensional  step  sizes  are 

Ar+  = 
Az+  = 
A  t+  = 


Ar 

CYL 

A £  _ 

L 

aA  t 

~w 


Nr 


M 


MAX 


(72) 

(73) 

(74) 


The  outer-radial  boundary  of  the  cylinder  is  always  at  non-dimensional  radius  of 
1.0.  Then,  for  the  two-material  problem,  non-dimensional  radius  values  in  the  sur¬ 
rounding  material  are  greater  than  1.0.  Similarly,  the  back  surface  of  the  cylinder 
is  always  at  non-dimensional  axial  dimension  of  1.0. 
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The  non-dimensional,  finite-difference  model  does  include  the  very  same 
non-dimensional  parameters  introduced  in  Section  2.2.  Three  additional  non- 
dimensional  parameters  in  the  finite-difference  model  are  the  non-dimensional  step 
sizes;  Ar,  Az  and  At.  One  can  think  of  the  index  variables  n,  m  and  j  as  being 
the  independent  variables  in  the  finite-difference  model,  while  the  non-dimensional 
step  sizes  are  additional  parameters  peculiar  to  the  finite-difference  model.  The 
continuous  model  should  be  equivalent  to  the  limiting  case  of  the  finite-difference 
model  as  the  step  sizes  approach  zero.  Two  other  property  ratios,  k/k  and  pcp/pc^, 
appear  in  the  non-dimensional  nodal  equations  for  nodes  on  the  boundary  between 
the  two  materials.  However,  these  two  property  ratios  are  functions  of  other  inde¬ 
pendent,  non-dimensional  parameters: 
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The  non-dimensional  nodal  equations  for  the  single-material  and  two-materi¬ 
al  problems  using  the  fully-insulated,  outer-radial  boundary  condition  are  given  in 
Appendices  A.l  and  A. 2,  respectively.  To  use  the  outer-radial  boundary  conditions 
for  maximum  heat  flux  out  or  maximum  heat  flux  in,  the  nodal  equations  for  nodes 
on  the  outer-radial  boundary  are  omitted,  and  the  temperatures  for  nodes  on  the 
outer-radial  boundary  are  set  to  the  appropriate  constant  values.  For  the  limiting 
case  of  maximum  heat  flux  out,  the  non-dimensional  temperature  at  the  outer- 
radial  boundary  is  zero.  For  the  limiting  case  of  maximum  heat  flux  in,  the  non- 
dimensional  temperature  at  the  outer-radial  boundary  is  zero  for  the  preheating 
problem  and  one  for  the  disturbance  problem  (see  Equations  27,  2S,  35  and  36). 


2.3.6  Stability  Criterion  When  using  an  explicit  method,  the  time  step  can 
not  be  chosen  arbitrarily.  The  time  step  must  be  small  enough  to  ensure  that  the 
system  is  stable.  The  instability  in  an  explicit  method  results  from  evaluating  the 


35 


finite-difference  expression  for  net  heat  flux  into  the  nodal  volume  at  the  old  time 
step  rather  them  at  the  new  time  step  or  a  combination  of  the  new  and  old  [4:p.215]. 
The  change  in  temperature  during  a  time  step  tends  to  reduce  the  temperature 
gradients  that  caused  the  net  heat  flux,  which  then  decreases  the  magnitude  of 
the  net  heat  flux  into  the  control  volume.  For  example,  if  heat  is  flowing  into  a 
control  volume  at  some  time  level,  the  temperature  in  that  control  volume  will 
increase.  The  increase  in  temperature  of  the  control  volume  will  tend  to  reduce 
the  temperature  gradients  that  caused  the  heat  to  flow  into  the  control  volume, 
which  in  turn  decreases  the  rate  at  which  heat  continues  to  flow  into  the  control 
volume. 

The  explicit  method  ignores  this  compensatory  effect.  Instead,  the  explicit 
method  evaluates  the  net  heat  flux  into  the  nodal  volume  at  the  old  time  level 
and  assumes  that  this  net  rate  of  energy  flow  into  the  control  volume  is  constant 
throughout  the  next  time  step.  It  then  predicts  the  temperature  at  the  node  for  the 
new  time  level  using  the  net  rate  of  energy  flow  evaluated  from  the  old  time  level. 
As  a  result,  if  the  time  step  is  too  large,  the  explicit  method  can  predict  a  change  in 
temperature  which  violates  the  second  law  of  thermodynamics.  The  overprediction 
in  each  step  of  the  explicit  method  can  also  cause  the  solution  to  oscillate  without 
necessarily  becoming  unstable.  Instabilities  occur  when  these  oscillations  grow  so 
that  eventually  the  temperatures  alternate  between  increasingly  larger  positive  and 
negative  values. 

Incropera  and  Dewitt  suggest  a  criterion  to  use  in  finding  the  maximum  al¬ 
lowable  time  step  [4 :p.215] .  The  criterion  they  suggest  is  easy  to  use  and  not  only 
assures  that  the  system  will  be  stable  but  also  helps  to  assure  that  the  solution  will 
not  oscillate.  For  heat  transfer  with  constant  boundary  conditions,  temperatures 
at  each  node  should  change  smoothly  in  one  direction.  Temperatures  shou'd  not 
increase  in  one  time  step,  decrease  in  the  next,  and  increase  in  the  following.  In¬ 
cropera  and  Dewitt  suggest  limiting  the  time  step  so  that  the  coefficients  in  front 
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of  the  temperature  terms  in  all  the  nodal  equations  (see  Appendix  A)  are  positive. 
It  is  easy  to  see  that  this  criterion  assures  stability;  if  this  criterion  is  satified, 
there  is  no  mechanism  for  non-dimensional  temperatures  to  become  negative.  It 
is  important  to  apply  the  criterion  to  all  nodal  equations  and  not  just  the  inte¬ 
rior  nodal  equations.  If  a  solution  diverges  at  any  node,  the  whole  solution  will 
eventually  diverge.  In  fact,  the  nodal  equations  for  nodes  on  the  boundaries  in  the 
two-dimensional  model  will  have  more  severe  stability  requirements  because  their 
dimensions  are  smaller.  Although  the  criterion  suggested  by  Incropera  and  Dewitt 
is  a  sufficient  but  not  necessary  criterion  for  stability,  the  time  steps  determined 
using  the  criterion  are  often  very  close  to  time  steps  which  cause  instability  in  the 
two-dimensional  model. 

For  the  single-material  problem,  the  nodal  equation  for  the  node  at  the  cen¬ 
terline  and  on  the  surface  with  the  larger  Biot  number,  usually  the  front  surface, 
has  the  most  severe  stability  requirement  (node  (a)  in  Figure  4): 

Ai+  <  (4  [i£7  N'YL  ~ 2M""  ”  (77) 

For  the  two-material  problem,  the  nodal  equation  with  the  most  severe  stability 
requirement  depends  on  the  values  of  the  property  ratios.  In  any  case,  the  node 
with  the  most  severe  stability  requirement  is  again  a  node  on  the  surface  with 
the  larger  Biot  number,  usually  the  front  surface.  The  node  with  the  governing 
requirement  may  be  the  node  at  the  centerline,  the  node  at  the  material  boundary 
or  any  of  the  nodes  between  the  material  boundary  and  the  outer-radial  boundary 
(node  (a),  node  (c)  or  nodes  in  region  (d)  in  Figure  5).  The  stability  requirement 
for  these  three  nodal  equations  in  their  respective  order  are 
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The  second  two  stability  requirements  can  become  more  severe  than  the  first  for 
cases  in  which  (pc?)  /  (pcj,)'  >  1  or  k'/k  >  1  or  both. 

2.8.7  Convergence  Criterion  for  the  Preheating  Problem.  The  preheating 
problem  models  the  time  after  the  instrumentation  is  turned  on  during  which 
the  Wheatstone  bridge  is  being  balanced.  The  preheating  of  the  quartz  cylinder 
continues  until  the  temperature  of  the  thin-film  sensor  stabilizes  enough  that  the 
bridge  can  be  balanced.  When  simulating  the  preheating  problem,  some  criterion 
must  be  used  to  determine  when  to  end  the  preheating  problem  and  begin  the 
disturbance  problem.  The  criterion  used  in  the  simulation  models  the  decision  to 
end  the  preheating  problem  by  quitting  when  the  rate  of  change  in  the  solution 
is  less  than  some  tolerance  value.  Specifically,  at  each  time  level  the  criterion 
picks  the  largest  change  in  non-dimensional  temperature  at  any  node,  divides  the 
absolute  value  of  that  largest  change  by  the  largest  non-dimensional  temperature 
in  the  system,  and  then  divides  this  normalized  maximum  change  by  the  time  step 
to  arrive  at  a  maximum  normalized  rate  of  change.  The  maximum  normalized 
rate  of  change  is  then  compared  with  the  tolerance  value.  Figure  8  illustrates  the 
decision. 

The  maximum  change  is  normalized  with  the  largest  non-dimensional  tem¬ 
perature  at  that  time  step  to  better  model  the  decision  process.  To  an  external 
observer  the  system  seems  to  stabilize  when  the  noticeable  rate  of  change  is  small. 
The  normalized  rate  of  change  measures  the  noticeable  rate  of  change.  For  exam¬ 
ple,  a  non-dimensional  temperature  change  of  0.002  in  one  time  step  can  represent 
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Figure  S.  Convergence  Criterion 


a  significant  rate  of  change  when  the  maximum  non-dimensional  temperature  value 
at  that  time  is  0.04.  However,  the  same  non- dimensioned  temperature  change  of 
0.002  in  another  time  step  can  represent  an  insignificant  and  unnoticeable  rate  of 
change  if  the  maximum  non-dimensional  temperature  value  at  that  time  is  on  the 
order  of  1. 


2.4  The  One- Dimensional  Series  Solution 


Thin-film  heat  transfer  gages  measure  the  temperature  on  the  front  surface 
of  the  gage  at  the  location  of  the  thin-film  sensor.  Surface  heat  flux  is  derived 
from  this  data  using  the  one-dimensional  series  solution  first  proposed  by  Cook 
and  Felderman  [2:pp. 561-562].  The  underlying  assumption  for  the  series  solution 
is  that  the  quartz  cylinder  behaves  as  an  initially  isothermal,  semi-infinite  solid. 


The  governing  equation  for  temperature  in  a  semi-infinite  solid  is  [4:p.202] 

(81) 


<w_  aw 

dt  dx 2 


where  x  is  the  distance  into  the  solid.  The  initial  condition  and  the  boundary 
condition  for  large  x  in  an  initially  isothermal,  semi-infinite  solid  are 


6(x,t  <  0)  =  0 
0(x  — *  00,  t  >  0)  =  0 


(82) 

(83) 


The  series  solution  is  based  on  the  solution  for  temperature  in  a  semi-infinite  solid 
with  constant  surface  temperature.  The  front  surface  boundary  condition  on  x  for 
the  problem  with  constant  surface  temperature  is 

6(x  =  0,t  >  0)  =  6,  (84) 

The  solution  for  temperature  in  the  semi-infinite  solid  with  constant  surface  tem¬ 
perature  is  [4:p.203] 

0(x,t)  =  6.c  rfc(^=)  (85) 
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Using  the  Duhammel  superposition  method  and  this  temperature  solution, 
one  can  derive  a  complete  expression  for  the  temperature  distribution  within  an 
initially  isothermal,  semi-infinite  solid  when  the  surface  temperature  is  some  known 
function  of  time.  This  expression  is  in  the  form  of  an  integration  over  time  of  the 
surface  temperature.  One  can  then  solve  for  the  surface  heat  flux  at  any  time  using 
Fourier’s  conduction  law  immediately  below  the  surface,  q,  =  —k^j  ^  .  The 
expression  for  surface  heat  flux  is  also  in  the  form  of  an  integration  over  time  of 
the  surface  temperature.  The  integral  expression  used  by  Cook  and  Felderman, 
modified  slightly  to  be  in  terms  of  basic  properties,  is 
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Cook  and  Felderman  proposed  to  approximate  the  exact  integral  solution  for 
surface  heat  flux  by  approximating  the  surface  temperature  function  as  a  piecewise 
linear  function  of  time.  Then,  the  integral  in  Equation  86  can  be  performed  exactly. 
The  resulting  equation  for  surface  heat  flux  is  in  the  form  of  a  series  summation 
of  surface  temperatures  at  discrete  values  of  time.  Rearranged  slightly,  Cook  and 
Felderman ’s  series  solution  for  surface  heat  flux  is  [5:p.l2] 
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For  a  solid  that  does  behave  like  an  initially  isothermal,  semi-infinite  solid, 
the  only  approximation  results  from  modeling  the  actual  surface  temperature  func¬ 
tion  as  a  piecewise  linear  function  of  time.  Cook  and  Felderman  tested  the  series 
solution  for  some  cases  in  which  surface  heat  flux  into  a  semi-infinite  solid  is  known 
analytically.  They  show  that  when  exact  values  for  surface  temperature  are  used  in 
the  series,  the  series  solution  is  well  behaved  and  fairly  accurate.  Surface  heat  flux 
values  tend  to  be  slightly  high,  but  accuracy  improves  when  smaller  time  intervals 
are  used. 
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2.5  Comparing  the  Series  Solution  Estimate  with  the  Simulated  Surface  Heat  Flux 

The  temperatures  from  the  two-dimensional  simulation  are  used  in  the  se¬ 
ries  solution  to  get  the  series  solution  estimate  for  external  surface  heat  flux.  The 
objective  of  this  analysis  is  to  evaluate  the  accuracy  of  the  series  solution  by  com¬ 
paring  these  series  solution  estimates  to  the  actual,  external  surface  heat  flux  from 
the  simulation. 


To  use  the  series  solution  with  temperatures  from  the  two-dimensional  model, 
the  series  solution  must  also  be  non-dimensionalized.  Using  the  non-dimensional 
form  for  the  disturbance  problem,  the  non-dimensional  series  solution  is 
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When  using  the  gages  in  the  laboratory,  one  calculates  a  single  value  for  the 
surface  temperature  from  the  change  in  resistance  of  the  thin-film  sensor.  When 
the  temperature  is  not  constant  over  the  front  surface  of  the  gage,  the  temperature 
determined  from  the  resistance  change  of  the  sensor  is  some  average  temperature 
of  the  surface  of  the  gage  at  the  thin  film. 


A  good  way  to  model  the  process  of  finding  a  single  average  temperature  is  to 
average  the  nodal  temperatures  in  the  region  covered  by  the  thin-film  sensor.  The 
average  should  be  a  weighted  average  in  which  the  individual  nodal  temperatures 
are  weighted  by  the  fraction  of  the  total  surface  area  of  the  sensor  covered  by 
each  of  the  nodal  regions.  For  the  two-dimensional  model,  a  simplified  form  of  the 
weighted  average  is 

a+  _  l/4fl\0  +  2E?°/**nfl+,,n 

(Ndisk  +  1  /2 ) 2 

The  average  surface  temperatures  can  then  be  used  in  the  series  solution  to  find 
the  series  solution  estimate  for  external  surface  heat  flux. 


The  surface  heat  flux  value  to  which  the  series  solution  estimate  is  compared 
is  the  actual  average,  external  surface  heat  flux  from  the  finite-difference  model 
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in  the  region  of  the  heated  disk.  From  Equation  (30),  the  non-dimensional,  ex¬ 
ternal  surface  heat  flux  at  each  node  in  the  region  of  the  heated  disk  during  the 
disturbance  problem  is 

1*i-  =  1  -  *+l.„  (80) 

The  average,  external  surface  heat  flux  is  found  using  Equation  (89)  and  replacing 
the  nodal  temperatures  with  the  nodal  surface  heat  fluxes. 

An  alternate  and  equivalent  method  for  finding  the  average,  external  surface 
heat  flux  from  the  finite-difference  model  is  to  use  the  average  surface  temperatures 
from  Equation  (89)  directly  in  Equation  (90).  Similarly,  an  alternate  and  equivalent 
method  for  finding  the  series  solution  estimate  is  to  use  the  temperatures  at  each 
node  in  the  non-dimensional  series  approximation  to  find  series  solution  estimates 
for  heat  flux  at  each  node  and  then  to  take  a  weighted  average  of  these  series 
solution  estimates. 

The  series  solution  and  finite-difference  estimates  for  external  surface  heat 
flux  are  compared  by  finding  their  difference.  The  difference  in  the  two  estimates  is 
automatically  normalized  with  respect  to  the  maximum  theoretical,  external  sur¬ 
face  heat  flux  by  the  non-dimensional  form  used  in  the  disturbance  problem.  The 
external  surface  heat  flux  is  maximum  if  the  surface  temperature  is  the  undisturbed 
initial  temperature,  T,: 

9*.  MAX  ~  front  ~  )  (01) 

This  maximum,  external  surface  heat  flux  is  1.0  in  non-dimensional  terms.  Then, 
the  difference  in  the  non-dimensional  series  solution  and  finite-difference  estimates 
multiplied  by  100.00  is  also  the  percent  difference  of  the  two  estimates  relative  the 
maximum  theoretical,  external  surface  heat  flux: 

Percent  Difference  =  {q+ER  —  <1+ D  )  x  100.00  (92) 

The  code  used  in  the  analysis  to  compare  the  series  solution  and  the  finite- 
difference  estimates  for  external  surface  heat  flux  is  given  in  Appendix  D.3. 
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III.  Check  Cases 


Check  cases  on  the  finite-difference  model  serve  the  same  purpose  that  cal¬ 
ibration  serves  in  laboratory  experiments.  Check  cases  test  the  model  for  gross 
errors  (i.e.  typing  errors  and  other  debugging  problems)  and  systematic  errors 
that  result  from  finite-differencing.  The  check  cases  will  each  be  simplified  prob¬ 
lems  for  which  an  analytical  solution  exists.  As  a  result,  they  will  check  individual 
aspects  of  the  model  separately.  The  check  cases  should  ultimately  show  how  ac¬ 
curate  the  model  is  for  its  intended  use  and  in  what  ways  or  under  what  conditions 
the  model  may  be  inaccurate. 

The  check  cases  for  the  finite-difference  model  are  divided  into  three  cate¬ 
gories.  The  first  two  are  steady-state  cases  for  one-dimensional  heat  transfer  in  the 
axial  direction  and  then  the  radial  direction.  The  third  catagory  includes  transient 
cases  for  one-dimensional  heat  transfer  in  the  axial  direction  only.  The  steady-state 
cases  primarily  check  the  accuracy  of  the  spatial  derivatives  in  each  of  the  two  di¬ 
rections.  The  accuracy  to  which  the  solutions  converge  on  the  correct  steady-state 
values  should  indicate  the  accuracy  of  the  derivatives.  Another  check  on  the  ac¬ 
curacy  of  the  derivatives  would  be  to  observe  the  transient  temperature  profiles 
during  the  steady-state  cases  to  assure  that  the  transient  solutions  progress  logi¬ 
cally  toward  the  steady-state  solution.  Also,  one  can  check  for  gross  errors  simply 
by  assuring  that  the  heat  transfer  is  indeed  one- dimensional  at  each  time  step. 

The  transient  cases  check  the  accuracy  of  the  time  derivatives  in  the  model 
in  conjunction  with  the  axial  derivatives  only.  However,  if  the  time  derivatives  are 
accurate  for  axial  heat  transfer,  then  they  should  also  be  accurate  for  radial  heat 
transfer.  The  transient  cases  also  check  the  accuracy  of  the  series  solution  for  sur¬ 
face  heat  flux  when  the  solid  does  behave  as  an  initially  isothermal,  semi-infinite 
solid.  Although  the  transient  cases  could  be  used  to  check  the  accuracy  of  tran¬ 
sient  temperature  profiles  in  the  solid,  the  transient  cases  were  used  to  check  the 
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accuracy  of  only  transient  surface  temperatures  since  the  model  is  primarily  used 
to  produce  transient  surface  temperatures  for  the  analysis.  The  series  solution  was 
checked  using  both  the  analytical  values  for  surface  temperature  and  the  surface 
temperatures  produced  by  the  finite-difference  model  in  order  to  distinguish  be¬ 
tween  inaccuracies  inherent  to  the  series  solution  and  those  caused  by  the  surface 
temperatures  produced  by  the  finite-difference  model. 

S.l  One- Dimensional,  Steady-State  Heat  Transfer  in  the  Axial  Direction 

The  one-dimensional,  steady-state  check  cases  for  heat  transfer  in  the  axial 
direction  are  separated  into  cases  using  the  preheating  problem  and  cases  using  the 
disturbance  problem.  To  obtain  one-dimensional  heat  transfer  in  the  preheating 
problem,  the  heat  generation  must  be  uniform  over  the  front  surface.  Uniform 
surface  heat  generation  is  obtained  by  making  the  heated  disk  as  wide  as  the 
cylinder  and  by  using  the  single-material  problem  only.  The  check  cases  for  the 
preheating  problem  include  cases  in  which  the  Biot  number  at  the  back  surface  is 
either  zero  or  non-zero. 

One-dimensional  heat  transfer  can  be  obtained  in  the  disturbance  problem 
using  either  uniform  surface  heat  generation  or  no  surface  heat  generation.  When 
surface  heat  generation  is  absent,  heat  transfer  will  be  one-dimensional  for  both 
the  single-material  problem  and  the  two-material  problem  with  property  ratios 
of  unity.  Because  the  results  for  the  two  problems  are  identical,  only  the  results 
for  the  single-material  problem  axe  shown.  The  check  cases  for  the  disturbance 
problem  include  cases  in  which  surface  heat  generation  is  either  uniform  over  the 
surface  of  the  cylinder  or  absent  and  in  which  the  Biot  number  at  the  back  surface 
is  either  zero  or  non-zero. 

The  analytical  solutions  for  one-dimensional,  steady-state  heat  transfer  in 
the  axial  direction  were  determined  using  thermal  networks  developed  from  the 
electrical  analogy  [4:pp. 64-65].  For  one-dimensional,  steady-state  heat  transfer  in 
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the  axial  direction,  Equation  (3)  reduces  to 


Heat  flux  is  constant  because  the  temperature  gradient  is  constant.  When  the  tem¬ 
perature  gradient  is  constant,  both  the  temperature  gradient  and  the  temperature 
profile  through  the  solid  are  completely  specified  by  the  temperatures  at  the  front 
and  back  surfaces.  Then,  for  one-dimensional,  steady-state  heat  transfer  in  the  ax¬ 
ial  direction,  the  following  system  of  simultaneous  equations  completely  specifies 
the  analytical  solution: 
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The  simultaneous  equations  governing  one- dimensional,  steady-state  heat 
transfer  in  the  axial  direction  are  similar  to  equations  governing  current  flow  in 
electrical  networks.  Therefore,  one  can  apply  the  mathematics  used  with  electric 
circuits  to  analyze  the  heat  transfer.  In  the  electical  analogy,  heat  flux  is  analogous 
to  current  density,  and  temperature  difference  is  analogous  to  voltage  difference. 
Nodes  in  the  thermal  network  represent  positions  whose  temperatures  appear  in 
the  simultaneous  equations.  The  four  nodes  in  the  thermal  networks  for  the  check 
cases  represent  the  external  fluid  at  the  front  and  back  surfaces  and  the  front  and 
back  surfaces  of  the  solid  itself.  Equation  (97)  is  a  statement  of  energy  balance  for 
the  front  and  back  surfaces  of  the  solid.  It  is  analogous  to  a  statement  of  Kirch- 
hoff’s  current  law  at  the  nodes  corresponding  to  the  front  and  back  surfaces  of  the 
solid.  Equations  (94)-(96)  are  thermal  analogies  to  Ohm’s  law  from  which  thermal 
resistances  for  convection  at  the  front  and  back  surfaces  and  conduction  through 


the  solid  can  be  defined  as  follows: 

(Rt  ,conv)  FRONT  ~  L  (98) 
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The  heat  generation  on  the  front  surface  is  modeled  as  a  constant  current  into  the 
network  at  the  front  surface. 

To  use  the  electrical  analogy  with  the  finite- difference  model,  non-dimensional 
thermal  resistances  and  surface  heat  generation  are  determined  by  non-dimension- 
alizing  the  terms  in  Equations  (94)-(97).  Non-dimensional  thermal  resistances, 
surface  heat  generation  and  fluid  temperatures  for  use  in  the  preheating  and  dis¬ 
turbance  problems  are  summarized  by  Table  1.  To  help  visualize  the  thermal 
networks,  the  initial  temperature,  T,,  can  be  defined  as  a  ground  potential  from 
which  all  other  potentials  are  referenced. 

Results  for  two  cases,  one  a  preheating  problem  and  the  other  a  disturbance 
problem  are  shown  in  Figures  9  and  10.  The  figures  include  a  diagram  of  the 
corresponding  thermal  network  for  the  conditions  of  the  check  case.  The  preheating 
problem  shown  in  Figure  9  uses  Biot  numbers  of  1.0  for  the  front  surface  and  0.5  for 
the  back  surface.  The  disturbance  problem  shown  in  Figure  10  uses  Biot  numbers 
of  0.5  for  the  front  surface  and  zero  for  the  back  surface  and  a  non-dimensional 
surface  heat  generation  of  0.2.  The  analytical  solution  for  this  case  is  independent 
of  the  Biot  number  on  the  front  surface.  All  cases  used  a  non-dimensional  time  step 
of  0.004  and  a  non-dimensional  step  size  of  0.1  in  the  axial  direction.  Results  for  the 
other  preheating  problems  and  disturbance  problems  are  given  in  Appendices  B.1.1 
and  B.1.2,  respectively. 

The  finite-difference  model  requires  a  tolerance  for  use  in  the  convergence 
criterion  for  all  the  steady-state  problems.  Section  2.3.7  explains  the  convergence 
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Table  1.  Summary  of  Non-Dimensional  Terms  in  the  Electrical  Analogy 


Disturbance  Problem: 


criterion  used  in  the  model.  The  finite-difference  model  converges  when  the  rate 
of  change  becomes  small  rather  than  when  the  error  in  the  steady-state  solution 
becomes  small.  As  a  result,  the  proximity  of  the  transient  solution  to  the  true,  an¬ 
alytical,  steady-state  solution  when  the  model  converges  varies  with  the  conditions 
as  well  as  with  the  tolerance  used  in  the  convergence  criterion.  The  steady-state 
check  cases  used  a  tolerance  of  0.01  in  the  convergence  criterion  (see  Section  2.3.7). 
For  all  check  cases  in  the  axial  direction,  a  tolerance  of  0.01  in  the  convergence  cri¬ 
terion  caused  the  transient  solution  to  converge  within  2  percent  of  the  steady-state 
solution. 

The  two  cases  shown  in  Figures  9  and  10  were  run  with  different  values  for  the 
tolerance  in  the  convergence  criterion  to  demonstrate  that  the  transient  solutions 
do  progress  logically  towards  the  steady-state  solution.  Figures  11  and  12  show  the 
results.  Also,  the  heat  transfer  in  all  cases  was  indeed  one-dimensional  indicating 
that  there  probably  are  no  gross  errors  in  the  derivatives  for  the  axial  direction. 

3.2  One- Dimensional,  Steady-State  Heat  Transfer  in  the  Radial  Direction 

There  axe  four  cases  for  one-dimensional,  steady-state  heat  transfer  in  the 
radial  direction.  To  obtain  one-dimensional  heat  transfer  in  the  radial  direction, 
the  heat  flux  across  the  front  and  back  surfaces  must  be  set  to  zero  by  setting 
the  Biot  numbers  and  the  surface  heat  generation  equal  to  zero.  Because  surface 
heat  flux  is  zero,  the  resulting  nodal  equations  in  both  the  preheating  and  the 
disturbance  problems  are  the  same. 

Both  the  single-material  and  the  two-material  problems  produce  one-dimen¬ 
sional  heat  transfer  in  the  radial  direction.  For  simplicity,  property  ratios  of  unity 
were  used  in  the  cases  run  on  the  two-material  problem.  Because  of  the  way 
in  which  the  radial  dimension  is  non-dimensionalized,  the  results  for  the  single¬ 
material  problem  and  the  two-material  problem  are  different.  Therefore,  results 
for  both  problems  are  given. 
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Preheating  Problem: 


Figure  li.  Sample  Result  for  Steady  State  Heat  Transfer  in  the  Axial  Direction 
Using  the  Preheating  Problem  with  Different  Values  for  the  Tolerance 


Disturbance  Problem 
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The  analytical  solutions  for  one-dimensional,  steady-state  heat  transfer  in 
the  radial  direction  were  determined  by  solving  the  differential  equation.  For  one¬ 
dimensional,  steady-state  heat  transfer  in  the  radial  direction,  Equation  (3)  in 
non-dimensional  form  reduces  to 


1  d  / 

r+  dr+  \  dr+  ) 


=  0 


Integrating  the  equation  twice  yields 

d0+  _  Ci 

5r+  r+ 
e+  =  C\  ln(r+)  +  Ci 


(101) 


(102) 

(103) 


where  C\  and  are  constants  whose  value  depends  on  the  boundary  conditions. 

When  the  inner-radial  boundary  condition  is  the  symmetry  condition  of 
Equation  (8),  C\  is  equal  to  zero.  When  the  inner-radial  boundary  condition 
is  a  specified  temperature,  it  is  necessary  in  general  to  specify  the  temperature  at 
an  inner  radius  other  than  zero,  so  one  or  more  of  the  inner-radial  nodes  must  be 
maintained  at  a  constant  temperature.  The  analytical  solution  for  the  rest  of  the 
model  is  then  equivalent  to  the  solution  for  heat  transfer  in  a  hollow  cylinder  with 
a  specified  temperature  at  the  inner  radius.  Heat  transfer  in  a  hollow  cylinder  with 
inner  radius 

=  l-5Ar+ 

is  simulated  by  equating  the  temperature  of  the  first  two  radial  nodes  to  a  con¬ 
stant  value.  Because  their  temperatures  are  held  constant,  the  nodal  equations  for 
these  first  two  radial  nodes  are  not  used.  Therefore,  the  check  cases  which  spec¬ 
ify  a  temperature  for  the  inner-radial  boundary  condition  do  not  check  the  radial 
derivatives  for  the  nodal  equations  on  the  centerline. 

The  first  radial  check  case  checks  the  radial  derivatives  for  nodes  on  the 
centerline  by  using  the  symmetry  condition  for  the  inner-radial  boundary  condition. 
The  next  three  radial  check  cases  check  the  three  limiting  cases  for  the  outer-radial 
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boundary  condition.  Table  2  summarizes  the  boundary  conditions  and  analytical 
solutions  for  the  four  one-dimensional,  steady-state  check  cases  for  radial  heat 
transfer. 

For  each  check  case  using  the  single-material  problem,  the  non-dimensional 
radial  step  size  was  0.1.  For  each  check  case  using  the  two-material  problem, 
the  non-dimensional  radial  step  size  was  0.2  and  the  non-dimensional  geometry 
parameter,  Rmax/R-cyl >  was  2-0.  All  cases  used  a  total  of  eleven  nodes  in  the 
radial  direction,  a  non-dimensional  time  step  of  0.003  and  a  tolerance  of  0.01. 
Table  3  summarizes  the  analytical  solutions  for  each  of  the  radial  check  cases  for 
both  the  single-material  and  two-material  problems. 

Figure  13  shows  the  results  for  the  second  check  case  using  the  single-material 
problem.  The  tolerance  used  in  the  convergence  criterion  was  again  varied  to 
demonstrate  that  the  transient  solution  does  progress  logically  toward  the  steady- 
state  solution,  and  this  result  is  shown  in  Figure  14.  The  results  for  the  other 
check  cases  using  the  single-material  and  two-material  problems  are  given  in  Ap¬ 
pendices  B.2.1  and  B.2.2,  respectively.  In  the  steady-state  cases  for  heat  transfer 
in  the  radial  direction,  some  of  the  solutions  converged  farther  from  the  true  steady 
stai,e  solution  than  others.  In  general,  cases  in  which  heat  flows  radially  outward 
seem  to  have  a  much  slower  rate  of  change  so  that  the  transient  solution  converges 
farther  from  the  true  steady-state  solution.  Cases  using  the  two-material  problem 
also  generally  converged  farther  from  the  true  steady-state  solution.  Nevertheless, 
the  transient  profiles  do  progress  logically  and  eventually  converge  on  the  steady 
state  solution.  The  heat  transfer  wras  indeed  one-dimensional  in  each  case  showing 
that  there  probably  are  no  gross  errors  in  the  radial  derivatives. 

3.3  One- Dimensional,  Transient  Heat  Transfer  in  the  Axial  Direction 

Three  cases  for  one-dimensional,  transient  heat  transfer  w’ere  used  to  check 
the  finite-difference  model  and  the  series  solution.  The  transient  check  cases  only 
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Boundary  Conditions  Analytical  Solution 
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Transfer  in  the  Radial  Direction 


Single-Material  Problem  Two-Material  Problem 
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Transfer  in  the  Radial  Direction 
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Figure  14.  Sample  Result  for  Steady  State  Heat  Transfer  in  the  Radial  Direction 
with  Different  Values  for  the  Tolerance 


include  problems  with  axial  heat  transfer.  This  should  be  sufficient  since  both  the 
axial  and  radial  derivatives  proved  accurate  for  steady-state  heat  transfer. 

The  transient  cases  were  used  in  the  following  ways: 

•  To  check  the  accuracy  of  the  finite-difference  model  in  producing  non-dimen¬ 
sional  surface  temperatures  and  surface  heat  flux  values. 

•  To  check  the  accuracy  of  the  non-dimensional  series  solution  for  cases  in 
which  the  solid  does  behave  as  an  initially  isothermal,  semi-infinite  solid. 

•  To  investigate  the  problems  encountered  when  using  the  non-dimensional 
surface  temperatures  produced  by  the  finite-difference  model  in  the  series 
solution. 

•  To  determine  suitable  spatial  and  time  steps  to  get  sufficient  accuracy  from 
the  model  in  the  analysis  and  to  estimate  this  accuracy. 

Since  heat  transfer  is  one- dimensional  in  these  check  cases,  it  is  not  necessary 
to  average  surface  temperatures  or  surface  heat  fluxes.  Nevertheless,  each  case 
averaged  the  temperatures  and  heat  fluxes  of  the  first  three  radial  nodes  on  the 
front  surface  in  order  to  check  the  equations  for  finding  average  surface  temperature 
and  heat  flux.  The  simulation  of  each  case  on  the  finite-difference  model  used  non- 
dimensional  spatial  step  sizes  of  0.1  and  three  different  non-dimensional  time  steps 
of  0.004,  0.001  and  0.0005. 

The  sections  that  follow  describe  each  of  the  transient  check  cases  giving 
the  analytical  solution  for  each  and  showing  how  each  was  implemented  on  the 
finite-difference  model.  The  final  section  summarizes  the  results  from  the  transient 
cases.  Graphs  from  the  first  transient  check  case— the  semi-infinite  solid  with 
convection — are  used  to  illustrate  points  in  the  discussion.  Graphs  for  the  other 
two  transient  check  cases  are  given  in  Appendices  B.3.1  and  B.3.2. 
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3.3.1  The  Semi- Infinite  Solid  with  Convection  One-dimensional  heat  trans¬ 
fer  in  a  semi-infinite  solid  with  constant  external  fluid  temperature  and  convection 
coefficient  is  readily  simulated  using  the  finite-difference  model  for  the  disturbance 
problem  with  a  fully-insulated,  outer-radial  boundary.  Figure  15  illustrates  the 
conditions  used  in  the  two-dimensional  model  to  simulate  this  transient  solution. 
The  parameters  7,  r  and  dBiBACK  are  set  to  zero.  Both  the  single-material  problem 
and  the  two-material  problem  work.  However,  since  both  produce  identical  results, 
the  results  from  only  one  are  shown. 

The  cylinder  will  behave  as  a  semi-infinite  solid  only  for  times  less  than  the 
time  it  takes  the  leading  edge  of  the  thermal  disturbance  to  travel  the  length  of 
the  cylinder.  Using  the  approximation  given  by  Equation  (2),  the  simulation  for 
this  transient  check  case  should  only  be  valid  for  non-dimensional  times  less  than 

t+ ~  4  =  0.0625  (104) 

lo 

The  check  case  was  run  slightly  longer  to  non-dimensional  time  of  0.075. 

The  analytical  solutions  for  non-dimensional  surface  temperature  and  surface 
heat  flux  for  this  transient  check  case  are  [4:pp.202-206] 

0+,  =  1  -  exp  ^(Bi)2t+"j  erfc  ((Bi)t+l^2^  (105) 

q+,  =  1  -  B+3  =  exp  ((5i)2f+)  erfc  ((Bf)t+1/2)  (106) 

where  Bi  is  the  Biot  number  at  the  front  surface  in  the  simulation,  which  was 
chosen  to  be  2.0. 

3.3.2  The  Semi-Infinite  Solid  with  Constant  Surface  Heat  Flux  One-dimen¬ 
sional  heat  transfer  in  a  semi-infinite  solid  with  constant  surface  heat  flux  is  simu¬ 
lated  by  a  modified  version  of  the  finite-difference  model  for  the  preheating  problem 
with  a  fully-insulated,  outer  radial  boundary.  By  setting  vBiFRONT  equal  to  zero 
and  the  ratio  RDISK / R-cyl  e0ual  to  one,  the  heat  flux  over  the  front  surface  of  the 


Semi-Infinite  Solid 


Two-Dimensional  Model 


Figure  15.  Simulating  a  Semi-Infinite  Solid  with  Constant  External  Fluid  Tem¬ 
perature  and  Convection  Coefficient 


62 


cylinder  is  constant  and  equal  to  the  surface  heat  generation: 

Qs  =  Qy  (107) 

Because  both  the  fluid  temperature  and  the  convection  coefficient  are  not  defined 
in  the  analytical  solution  for  this  case,  the  non-dimensional  form  of  the  dependent 
variables  must  be  changed  to 


0++  _  T~Ti 

P  ( q,L/k ) 

(108) 

p9++  =  — 

9a 

(109) 

The  non-dimensional  equations  using  this  non-dimensional  form  for  the  dependent 
variables  are  again  identical  to  the  previous  equations  except  in  the  boundary 
conditions  for  the  front  and  back  surfaces.  The  nodal  equations  for  nodes  on  the 
front  and  back  surfaces  can  be  easily  modified  to  use  this  non-dimensional  form  in 
the  finite-difference  model. 

Figure  1G  illustrates  the  conditions  used  in  the  two-dimensional  model  to 
simulate  this  transient  solution.  The  Biot  number  for  the  back  surface  is  set  to 
zero.  Because  the  simulation  uses  surface  heat  generation,  only  the  single-material 
problem  will  produce  one-dimensional  heat  transfer.  The  simulation  is  again  only 
valid  for  times  less  than  the  time  it  takes  the  leading  edge  of  the  disturbance  to 
travel  the  length  of  the  cylinder.  This  transient  check  case  was  also  run  to  non- 
dimensional  time  of  0.975. 


The  analytical  solutions  for  non-dimensional  surface  temperature  and  surface 
heat  flux  using  the  modified,  non-dimensional  form  for  the  preheating  problem  are 
[4:pp. 202-206] 


*.++  = 


2 1+]/2 


1/2 


,++  - 


=  1 


(HO) 

(111) 


Because  the  surface  heat  flux  is  inherent  in  the  finite-difference  model,  it  is  not 
necessary  to  check  the  finite-difference  model’s  estimate  for  surface  heat  flux.  This 
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Semi-Infinite  Solid 


Two-Dimensional  Model 
(Single-Material  Problem  Only) 
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transient  check  case  is  still  useful,  however,  for  checking  the  accuracy  of  the  finite- 
difference  model’s  estimate  for  surface  temperature  and  the  accuracy  of  the  series 
solution’s  estimates  for  surface  heat  flux  using  both  the  analytical  surface  temper¬ 
atures  and  the  surface  temperatures  produced  by  the  model. 

S.S.S  The  Plane  Wall  with  Convection  One-dimensional  heat  transfer  in 
a  plane  wall  with  constant  external  fluid  temperature  and  convection  coefficient 
on  both  sides  is  also  easily  simulated  using  the  disturbance  problem  with  a  fully- 
insulated,  outer-radial  boundary  condition.  The  parameter  7  is  set  to  zero,  while 
the  parameter  r  is  set  to  one.  The  Biot  numbers  for  the  front  and  back  surfaces 
must  be  equal.  Because  surface  heat  generation  is  not  used,  both  the  single¬ 
material  and  two-material  problems  work  for  this  transient  check.  However,  the 
results  are  again  identical,  so  the  results  from  only  one  are  shown.  Figure  17  illus¬ 
trates  the  conditions  used  in  the  model  to  simulate  one-dimensional  heat  transfer 
in  a  plane  wall  with  convection. 

Although  the  simulation  for  the  heat  transfer  in  a  plane  wall  is  valid  for  all 
time,  the  cylinder  ceases  to  behave  as  a  semi-infinite  solid  for  times  greater  than 
the  time  it  takes  the  leading  edge  of  the  disturbance  to  reach  the  midplane  of  the 
wall: 

t+  %  1/64  «  0.016  (112) 

For  times  greatei  than  this,  then,  one  can  not  use  this  transient  check  case  to 
check  the  accuracy  of  the  series  solution.  However,  this  transient  check  case  is  still 
useful  for  checking  the  accuracy  of  values  for  surface  temperature  and  heat  flux 
produced  by  the  finite-difference  model.  The  transient  check  case  using  the  plane 
wall  solution  was  also  run  to  non-dimensional  time  of  0.075. 

The  analytical  solutions  for  the  plane  wall  are  infinite-series  solutions.  The 
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Figure  If.  Simulating  the  Plane  Wall  with  Constant  External  Fluid  Temperature 
and  Convection  Coefficient  on  Both  Sides 


66 


solutions  can  be  approximated  by  the  the  first  four  terms  to  give  [4:pp.  183-187] 
Of  =  1  -  ^C'/exp  (-4C/2<+)  cos(O)  (113) 

/=i 

qf  =  l-^J+  =  X:C’iexp(-4Ci2t+)cos(0)  (114) 


i=i 


where  the  coefficient  Ct  is  given  by 

Ci  = 


4  sin  (0) 


20  + sin  (20) 

and  the  constants  Q  are  the  Ith  eigenvalues  of  the  characteristic  equation 

2 0  tan  (0)  =  Bi 


(115) 


(116) 


Bi  is  the  Biot  number  for  both  front  and  back  surfaces.  For  the  Biot  number  of 
2.0  used  in  the  simulation,  the  first  four  eigenvalues  are 

Ci 
C2 
C3 

(4 

For  non-dimensional  time  of  0.002,  the  magnitude  of  the  fourth  term  in  the  series 
is  0.0104.  Since  this  is  an  alternating  series,  the  error  in  the  solution  using  the 
truncated  series  is  less  than  0.0104  in  magnitude  for  all  non-dimensional  times 
greater  than  0.002. 


0.8603 

rad 

(117) 

3.4256 

rad 

(118) 

6.4376 

rad 

(119) 

9.5293 

rad 

(120) 

S.S.4  Results  from,  the  Transient  Check  Cases  Because  the  first  check  case 
is  similar  to  the  disturbance  problem  of  the  analysis,  the  model  should  behave 
similarly  in  the  analysis  as  in  this  first  check  case.  For  this  reason,  the  accuracy 
of  the  model  in  the  first  check  case  should  give  a  good  indication  of  the  accuracy 
of  the  model  in  the  analysis.  The  one  major  shortcoming  is  the  fact  that  this  first 
check  case  does  not  exercise  the  radial  derivatives.  Nevertheless,  the  results  from 
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the  first  transient  check  case  are  given  in  this  section,  and  from  these  results  an 
estimate  of  the  expected  accuracy  of  the  model  in  the  analysis  is  determined.  The 
graphs  of  the  results  for  the  second  and  third  transient  check  cases  are  given  in 
Appendices  B.3.1  and  B.3.2,  respectively.  For  the  most  part,  the  trends  in  the 
accuracy  seen  in  the  first  check  case  are  also  seen  in  the  other  two  transient  check 
cases. 

Figure  18  shows  the  surface  temperatures  produced  by  the  finite-difference 
model  in  the  first  check  case  using  values  of  0.004,  0.001,  and  0.005  for  the  non- 
dimensional  time  step.  There  are  some  problems  with  the  surface  temperatures 
produced  by  the  finite-difference  model.  First  of  all,  as  the  non-dimensional  time 
step  is  reduced,  the  surface  temperatures  converge  on  values  lower  than  the  analyt¬ 
ical  values.  This  error  in  the  surface  temperatures  produced  by  the  finite-difference 
model  is  most  likely  related  to  the  magnitude  of  the  temperature  gradients  near 
the  surface.  The  error  is  greatest  for  small  times  when  the  temperature  gradi¬ 
ents  are  largest.  Nevertheless,  the  surface  temperatures  produced  by  the  model 
are  quite  accurate  if  sufficiently  small  non-dimensional  spatial  and  time  steps  are 
used.  Using  spatial  steps  of  0.1  and  a  time  step  of  0.001  or  less,  the  finite-difference 
estimates  for  surface  temperature  are  within  2.7  percent  of  the  analytical  value  at 
non-dimensional  time  of  0.01  and  within  0.03  percent  of  the  analytical  value  at 
non-dimensional  time  of  0.05. 

The  rate  of  change  of  surface  temperature  with  respect  to  time  as  seen  by  the 
slope  of  the  curves  in  Figure  18  is  severely  affected  by  these  seemingly  small  errors 
in  the  surface  temperatures  produced  by  the  model.  For  small  time,  the  slope  of 
the  curves  for  the  surface  temperatures  produced  by  the  model  is  significantly  lower 
than  in  the  analytical  solution.  For  larger  times,  the  slope  is  slightly  higher.  This 
error  in  the  rate  of  change  of  surface  temperature  is  important  because  the  series 
solution  uses  the  rate  of  change  of  surface  temperature  rather  than  the  surface 
temperatures  themselves  to  determine  the  surface  heat  flux. 
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Figure  IS. 


Finite-Difference  Model  Values  for  df  vs.  t+  for  the  Semi-Infinite  Solid 
with  Convection 


Stability  and  oscillation  are  also  a  concern  with  the  explicit  finite-difference 
model.  Oscillations  can  have  a  significant  effect  on  the  results  for  surface  heat  flux 
using  the  series  solution  because  the  oscillations  cause  large  errors  in  the  calculated 
rate  of  change  of  surface  temperature.  Figure  18  shows  oscillations  in  the  finite- 
difference  estimates  for  surface  temperature  using  a  non-dimensional  time  step  of 
0.004. 

Oscillations  are  a  more  likely  problem  for  small  time  when  the  curvature  of 
the  surface  temperature  function  of  time  is  greatest.  The  oscillations  are  caused  by 
the  over  prediction  in  the  explicit  scheme  as  described  in  Section  2.3.6.  The  non- 
dimensional  surface  temperature  at  the  end  of  the  first  time  step  is  grossly  over 
predicted  because  the  non-dimensional  time  step  is  too  big  to  accurately  follow 
the  solution  when  the  rate  of  change  in  slope  is  high.  The  over  predicted  surface 
temperature  at  the  end  of  the  first  time  step  causes  the  rate  of  change  for  the  next 
time  step  to  be  low,  which  in  turn  causes  the  surface  temperature  at  the  end  of 
the  following  time  step  to  be  under  predicted.  The  process  repeats  itself  a  few 
times  until  the  oscillations  die  out  when  the  curvature  of  the  surface  temperature 
function  becomes  small  enough  for  the  finite-difference  model  to  follow  the  solution 
with  the  given  time  step.  As  is  readily  seen  in  Figure  18,  the  oscillations  cause 
much  greater  error  in  the  slope  of  the  curve  for  surface  temperature  than  in  the 
actual  values  for  surface  temperature.  The  time  step  needed  to  prevent  oscillations 
depends  on  the  conditions  in  the  simulation  and  the  spatial  step  sizes  being  used. 

Because  of  the  way  in  which  the  disturbance  problem  is  non-dimensionalized, 
the  accuracy  of  the  finite-difference  estimates  for  surface  heat  flux  depends  solely  on 
the  accuracy  of  the  surface  temperatures  produced  by  the  model  (see  Equation  (90) 
and  the  discussion  in  Section  2.5).  Figure  19  shows  the  finite-difference  estimates 
for  surface  heat  flux  from  the  first  transient  check  case  using  non-dimensional 
time  steps  of  0.004,  0.001,  and  0.0005.  The  finite-difference  estimates  tend  to  be 
high  since  the  surface  temperatures  tend  to  be  low.  Also,  if  the  non-dimensional 


70 


time  step  is  too  big,  the  finite-difference  estimates  for  surface  heat  flux  oscillate 
about  the  analytical  values.  Nevertheless,  the  accuracy  is  quite  good  when  using 
small  enough  step  sizes.  Using  non-dimensional  spatial  step  sizes  of  0.1  and  a  non- 
dimensional  time  step  of  0.001  or  less,  the  finite-difference  estimate  for  surface  heat 
flux  is  within  0.7  percent  of  the  analytical  value  at  non-dimensional  time  of  0.01 
and  within  0.02  percent  of  the  analytical  value  at  non-dimensional  time  of  0.05. 

Errors  in  the  series  solution  for  surface  heat  flux  when  using  surface  tempera¬ 
tures  produced  by  the  finite-difference  model  may  be  due  either  to  the  inaccuracies 
in  the  surface  temperatures  produced  by  the  model  or  to  the  inherent  inaccuracy 
of  the  series  solution.  To  distinguish  between  the  two  causes,  the  series  solution 
for  surface  heat  flux  was  run  using  both  the  analytical  values  for  surface  temper¬ 
atures  and  the  surface  temperatures  produced  by  the  finite-difference  model.  It 
is  important  to  ensure  that  large  errors  in  the  series  solution  during  the  analysis 
are  not  caused  by  inaccuracies  in  the  surface  temperatures  produced  by  the  finite- 
difference  model.  If  this  is  the  case,  then  the  finite-difference  model  does  not  model 
the  laboratory  experiment  with  sufficient  accuracy,  and  the  results  are  meaning¬ 
less.  Running  the  series  solution  with  the  analytical  values  for  surface  temperature 
is  also  helpful  for  evaluating  any  inherent  inaccuracy  in  the  series  solution.  The 
series  solution  is  not  an  exact  solution,  and  its  accuracy  depends  on  the  time  step 
used  in  the  series. 

Figure  20  shows  the  series  solution  estimates  for  surface  heat  flux  in  the  first 
transient  check  case  when  using  analytical  values  for  surface  temperature.  Because 
the  approximation  in  the  series  solution  comes  from  assuming  a  piecewise  linear 
function  of  time  for  the  surface  temperature,  the  series  solution  is  more  prone 
to  error  at  small  time  when  the  curvature  in  the  surface  temperature  function  of 
time  is  larger.  Since  more  recent  terms  in  the  series  solution  are  weighted  more 
heavily,  the  series  estimate  becomes  increasingly  more  accurate  for  larger  time  as 
the  curvature  in  the  surface  temperature  function  of  time  in  this  transient  check 
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Scries  Estimate  for  qf 
(Using  Analytical  Values  for  6 * ) 


vs. 


t* 


Figure  20.  Series  Solution  Estimates  for  qf  (Using  Analytical  Values  for  9f)  vs. 
t+  for  the  Semi-Infinite  Solid  with  Convection 
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case  decreases.  For  cases  in  which  the  surface  temperature  rises  monotonically  such 
as  this,  the  series  solution  will  tend  to  overestimate  surface  heat  flux.  Overall,  the 
series  solution  is  quite  accurate  if  a  small  enough  time  step  is  used  in  the  solution. 
Using  a  non-dimensional  time  step  of  0.001  or  less  in  the  series  solution  gave  series 
solution  estimates  for  surface  heat  flux  within  0.6  percent  of  the  analytical  value 
at  non-dimensional  time  of  0.01  and  within  0.05  percent  of  the  analytical  value  at 
non-dimensional  time  of  0.05. 

There  are,  however,  significant  problems  with  the  series  solution  when  us¬ 
ing  surface  temperatures  produced  by  the  finite-difference  model.  These  errors 
are  caused  by  inaccuracies  in  the  surface  temperatures  and  not  by  the  series  so¬ 
lution.  The  error  in  the  series  solution  estimate  is  significant  even  though  the 
error  in  the  value  of  the  surface  temperatures  is  small  because  the  series  solution 
uses  the  rate  of  change  of  surface  temperature  to  evaluate  surface  heat  flux.  In 
comparison,  the  finite- difference  estimates  for  surface  heat  flux  do  not  have  the 
same  problems  because  the  finite-difference  estimates  use  the  actual  value  of  the 
surface  temperature  rather  than  the  rate  of  change  of  surface  temperature.  First 
of  all,  the  rate  of  change  of  surface  temperatures  produced  by  the  finite-difference 
model  is  low  for  small  time,  so  the  series  estimates  are  also  initially  low.  For  larger 
times,  the  rate  of  change  is  slightly  high  which  may  cause  the  series  solution  to 
be  higher  than  it  would  be  otherwise.  However,  any  over  prediction  in  the  series 
solution  for  larger  times  may  be  offset  by  the  initial  under  prediction  as  is  the  case 
in  the  second  transient  check  case  (see  Figur:  57  in  Appendix  B.3.1).  Secondly, 
if  the  non-dimensional  time  step  used  in  the  model  is  too  big  so  that  oscillations 
are  present  in  the  surface  temperatures,  the  series  solution  estimates  can  oscillate 
severely  due  to  the  extreme  oscillation  in  the  rate  of  change  of  surface  temperature. 

Figure  21  shows  the  series  solution  estimates  for  surface  heat  flux  when  us¬ 
ing  the  surface  temperatures  produced  by  the  finite-difference  model  in  the  first 
transient  check  case  with  non-dimensional  time  steps  of  0.004,  0.001  and  0.0005  in 
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both  the  finite-differnce  model  and  the  series.  Even  with  a  non-dimensional  time 
step  of  0.0005,  the  initial  under  prediction  in  the  series  solution  is  severe.  A  severe 
initial  under  prediction  occurs  in  all  three  transient  check  cases  (see  Figures  57 
and  61  in  Appendices  B.3.1  and  B.3.2,  respectively).  However,  by  non-dimensional 
time  of  0.01,  the  series  solution  regains  its  accuracy.  Using  non-dimensional  spatial 
step  sizes  of  0.1  in  the  finite-difference  model  and  a  non-dimensional  time  step  of 
0.001  or  less  in  both  the  finite-difference  model  and  the  series  solution  gave  series 
solution  estimates  for  surface  heat  flux  within  2.7  percent  of  the  analytical  value 
at  non-dimensional  time  of  0.01  and  within  0.7  percent  of  the  analytical  value  at 
non-dimensional  time  of  0.05. 

In  summary,  the  finite-difference  model  produces  surface  temperatures  with 
good  accuracy  usin~  non- dimensional,  spatial  step  sizes  less  than  or  equal  to  0.1 
and  a  non-dimensional  time  step  less  than  or  equal  to  0.001.  The  series  solution 
is  also  quite  accurate  given  exact  data  when  the  non-dimensional  time  step  in 
the  series  is  less  than  or  equal  to  0.001.  However,  the  series  solution  for  surface 
heat  fiux  using  the  surface  temperatures  produced  by  the  finite-diffe.nce  model 
does  have  some  problems.  For  small  times,  the  series  solution  under  predicts  the 
surface  heat  flux  due  to  the  large  error  in  the  initial  rate  of  change  of  the  surface 
temperatures  produced  by  the  finite-difference  model.  The  series  solution  can 
also  have  significant  errors  if  oscillations  are  present  in  the  surface  temperatures 
produced  by  the  finite-difference  model.  Nevertheless,  the  series  solution  for  surface 
heat  flux  using  surface  temperatures  produced  by  the  model  is  accurate  for  non- 
dimensional  times  greater  than  0.01  if  non-dimensional  spatial  step  sizes  less  than 
or  equal  to  0.01  are  used  in  the  model  and  non-dimensional  time  steps  less  than  or 
equal  to  0.001  are  used  in  both  the  model  and  the  series. 

Accurate  results  in  the  analysis  require  accurate  results  for  both  the  finite- 
difference  estimates  for  surface  heat  flux  and  the  series  solution  estimates  for  surface 
heat  flux  using  the  surface  temperatures  produced  by  the  finite-difference  model. 
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For  non-dimensional  times  greater  than  0.01  in  the  first  transient  check  case,  the 
finite-difference  estimates  are  within  0.7  percent  of  the  analytical  values  for  surface 
heat  flux  and  the  series  solution  estimates  are  within  2.7  percent  of  the  analytical 
values  when  using  non-dimensional  spatial  step  sizes  of  0.1  in  the  finite-difference 
model  and  non-dimensional  time  steps  less  than  or  equal  to  0.001  in  both  the  model 
and  the  series.  To  be  conservative,  one  can  assume  that  these  percentages  apply  to 
the  maximum  non-dimensional,  external  surface  heat  flux  of  1.0.  Then,  the  non- 
dimensional  error  in  the  finite-difference  estimate  for  surface  heat  flux  is  always  less 
than  0  007  and  the  non-dimensional  error  in  the  series  solution  estimate  is  always 
less  than  0.027.  This  error  estimate  is  more  conservative  for  longer  times  since  the 
error  decreases  with  time.  It  is  important  to  note  that  the  transient  check  cases 
did  not  employ  the  radial  derivatives.  One  would  expect  the  radial  derivatives  to 
add  some  additional  error  to  the  results  in  the  analysis.  Then,  a.  reasonable  and 
conservative  estimate  of  the  expected  error  in  the  percent  difference  used  in  the 
analysis,  (?ser  —  <?f  d  ) x  100  00  •  for  non-dimensional  times  greater  than  0.01  would 
be  4.0.  Figure  22  shows  the  finite-difference  and  the  series  solution  estimates  for 
surface  heat  flux  using  the  surface  temperatures  produced  by  the  finite-difference 
model  in  the  first  transient  check  case  with  a  non-dimensionai  time  step  of  0.0005 
in  both  the  model  and  fhe  series. 

The  transient  check  cases  also  point  out  the  need  to  collect  temperature  data 
carefully  when  using  the  series  solution  on  experimental  data.  The  series  solution 
is  extremely  sensitive  to  enors  in  the  data  for  surface  temperatures  because  the 
series  solution  uses  the  rate  of  change  of  temperature  to  evaluate  surface  heat  flux. 
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IV.  Analysis  and  Results 


This  analysis  investigates  the  effect  of  various  parameters  in  the  two-dimen¬ 
sional  model  on  the  accuracy  of  the  series  solution  for  the  average  surface  heat 
flux  in  the  region  of  the  heated  disk.  Because  the  total  number  of  parameters  in 
the  model  is  large,  the  analysis  is  limited  by  allowing  only  some  of  the  parameters 
to  vary  and  setting  the  others  constant.  A  limited  range  or  number  of  values  is 
investigated  for  each  of  the  parameters  which  is  allowed  to  vary. 

Section  2.2  lists  the  non-dimensional  parameters  in  the  two-dimensional  mod¬ 
el.  Four  additional  non-dimensional  parameters  occur  in  the  finite-difference  model. 
They  are  the  tolerance  used  in  the  convergence  criterion  for  the  preheating  problem 
and  the  three  non-dimensional  step  sizes;  Az+,  Ar+  and  Af+. 

Some  of  the  parameters  are  set  constant  for  all  runs.  The  Biot  numbers  at  the 
front  and  back  surfaces  during  the  preheating  problem  are  set  to  0.01.  This  value 
is  chosen  using  a  typical  convection  coefficient  for  free  convection  of  4.0  W/m2  •  K 
and  typical  values  for  L  and  k  from  References  [3]  and  [7:p.672],  respectively.  The 
parameter  r  in  the  disturbance  problem  is  set  to  zero,  and  the  Biot  number  at 
the  back  surface  during  the  disturbance  problem  is  kept  at  0.01.  Setting  r  to 
zero  and  keeping  dBiBACK  equal  to  0.01  implies  that  the  conditions  at  the  back  of 
the  test  specimen  remain  the  same  for  both  the  preheating  and  the  disturbance 
problems  and  are  typical  conditions  for  free  convection.  Only  two  values,  0.1  and 
1.0,  are  used  to  investigate  a  reasonable  range  for  the  Biot  number  at  the  front 
surface  during  the  disturbance  problem.  Because  pBiFRONT  is  held  constant,  the 
parameter  /?  equivalently  describes  the  Biot  number  at  the  front  surface  during  the 
disturbance  problem  when  preheating  is  used. 

The  transient  check  cases  show  that  the  percent  difference,  {q+ER  —  qFD  )  x 
100.00  ,  should  be  accurate  to  within  4.0  in  the  analysis  for  non-dimensional 
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times  greater  than  0.01  when  using  non-dimensional  spatial  step  sizes  of  0.1  or  less 
and  a  non-dimensional  time  step  of  0.001  or  less.  The  spatial  step  sizes  actually 
used  depend  on  the  value  of  the  geometry  parameter  Rdisk/RCyl  •  The  geometry 
parameter  Rdisk/RCyl  ^  to  1/6,  1/8  and  1/12  using  non-dimensional  spatial 
step  sizes  of  0.111  (1/9),  0.083  (1/12)  and  0.056  (1/18),  respectively.  The  largest 
non-dimensional  time  step  used  in  any  of  the  runs  is  0.001. 

The  time  step  required  to  maintain  stability  in  the  finite-difference  model 
varies  between  runs.  In  some  cases,  a  time  step  much  smaller  than  0.001  is  used. 
The  accuracy  of  the  series  solution  improves  when  smaller  time  steps  are  used  in 
the  series.  However,  each  run  uses  the  temperature  values  at  non-dimensional  time 
steps  of  0.001  in  the  series  solution  regardless  of  the  time  step  used  in  the  finite- 
difference  model  to  avoid  disturbing  the  results  by  varying  the  time  step  used  in  the 
series  solution.  All  runs  are  ended  at  non-dimensional  time  of  0.07  since  thin-film 
gages  are  generally  not  used  for  times  longer  than  the  estimated  time  it  takes  the 
leading  edge  of  a  thermal  disturbance  to  reach  the  back  surface  of  the  cylinder, 
which  is  approximately  0.0625  in  non-dimensional  time  (see  Equation  (2)). 

In  another  attempt  to  clarify  and  simplify  the  analysis,  the  runs  are  separated 
into  groups  which  isolate  the  causes  of  two-dimensional  heat  transfer.  One  group 
of  runs,  the  adiabatic  cases,  investigates  the  effect  of  localized  heat  generation  and 
preheating  with  no  heat  flux  across  the  outer-radial  boundary  of  the  cylinder.  An¬ 
other  group  of  runs,  the  non-adiabatic  cases  with  no  heat  generation  or  preheating, 
investigates  the  effect  of  allowing  heat  flux  across  the  outer-radial  boundary  of  the 
cylinder  in  the  absence  of  heat  generation  or  preheating.  A  third  group  of  runs, 
the  combined  cases  or  non-adiabatic  cases  with  heat  generation  and  preheating,  in¬ 
vestigates  the  effect  of  adding  preheating  and  subsequent  heat  generation  to  a  few 
of  the  runs  from  the  non-adiabatic  cases. 
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4-1  Adiabatic  Cases 

The  adiabatic  cases  investigate  the  effect  on  the  series  solution  of  the  two- 
dimensional  heat  transfer  caused  by  heat  generation  and  preheating  in  the  absence 
of  heat  flux  across  the  outer-radial  boundary  of  the  cylinder.  All  runs  in  the 
adiabatic  cases  use  the  single  material  problem  with  the  fully-insulated,  outer-radial 
boundary  condition  cf  Equations  (26)  and  (34).  All  runs  also  use  the  preheating 
problem  except  the  few  which  investigate  the  effect  of  heat  generation  during  the 
disturbance  with  no  prior  preheating. 

A  number  of  parameters  are  allowed  to  vary  for  the  adiabatic  cases.  Values 
between  0.0  and  0.5  are  investigated  for  the  parameter  7.  The  geometry  parameter 
L/Rcyl  is  varied  between  0.2  and  5.0,  while  the  geometry  parameter  Rd1sk/RCYl 
is  varied  between  the  three  values  1/6,  1/8  and  1/12.  The  tolerance  used  in  the 
convergence  criterion  for  the  preheating  problem  is  varied  between  the  three  values 
of  0.5,  0.1  and  0.05;  and  a  few  runs  do  not  use  any  preheating.  Finally,  all  runs 
are  duplicated  using  values  of  10.0  and  100.0  for  the  parameter  /?. 

The  results  from  the  adiabatic  cases  show  that  preheating  and  subsequent 
surface  heat  generation  should  have  no  effect  on  the  results  as  long  as  the  transients 
from  the  preheating  problem  are  allowed  to  settle  down  before  beginning  the  test. 
Figure  23  shows  the  results  with  the  parameter  7  varying  and  using  a  value  of  10.0 
for  /?,  1/8  for  RDlSK / HCKi,  1-0  for  L/Rcyl  and  0.1  for  the  tolerance.  The  results 
seem  to  show  a  trend  in  the  series  solution  to  overestimate  slightly  with  increasing 
value  of  7.  However,  the  percent  differences  are  less  than  the  expected  error  in 
the  model.  Using  a  value  of  0.5  for  7,  the  percent  difference  is  between  0.5  and 
0.7  between  non-dimensional  time  of  0.01  and  0.07,  and  the  percent  difference  is 
less  using  smaller  values  for  7.  When  a  larger  value  for  (3  is  used,  the  percent 
difference  is  larger  but  still  less  than  the  expected  error  in  the  model.  Figure  63  in 
Appendix  C.l  shows  the  results  with  7  varying  using  a  value  of  100.0  for  /3.  The 
percent  difference  using  a  value  of  0.5  for  7  increases  to  between  2.0  and  3.0  when 
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Preheating 


0  is  100.0. 

Although  heat  generation  does  not  affect  the  results  when  the  transients  from 
the  preheating  problem  axe  allowed  to  settle  down,  the  effect  of  heat  generation  is 
large  if  the  preheating  problem  is  not  run.  Figure  24  shows  the  results  with  the 
tolerance  for  the  convergence  criterion  of  the  preheating  problem  varying  and  using 
a  value  of  10.0  for  0,  1/8  for  RD1Sh-l Rcyli  1-0  f°r  L/Rcyl  and  values  of  0.05  and 
0.2  for  7.  Tolerances  of  0.5,  0.1  and  0.005  cause  the  preheating  problem  co  con¬ 
verge  at  non-dimensional  times  of  0.118,  2.600  and  10.326,  respectively.  Changing 
the  tolerance  between  these  three  values  cause  only  slight  changes  in  the  percent 
differences  which  remained  below  the  expected  error  of  the  model.  However,  there 
is  significant  overestimation  when  heat  generation  is  used  with  no  preheating.  The 
percent  difference  when  using  a  value  of  0.2  for  7  and  without  preheating  is  be¬ 
tween  6.0  and  10.0.  Figure  64  in  Appendix  C.l  shows  even  greater  overestimation 
when  using  a  value  of  100.0  for  0. 

The  overestimation  in  the  series  solution  caused  by  not  allowing  the  transients 
to  settle  down  during  the  preheating  problem  stands  to  reason.  The  overestima¬ 
tion  results  from  greater  increases  in  surface  temperature  than  would  be  present 
from  the  external  disturbance  alone  since  the  disturbance  from  the  heat  generation 
has  not  established  itself  yet.  The  series  solution  attributes  the  greater  increase  in 
surface  temperature  to  a  greater  external  surface  heat  flux  than  is  actually  present. 
It  is  likely,  however,  that  if  one  used  a  thin-film  gage  for  very  short  times,  i.e.  only 
as  long  as  the  assumption  of  one-dimensional  heat  transfer  was  valid,  and  were 
able  to  expose  the  system  to  the  disturbance  at  the  same  time  as  turning  on  the 
instrumentation,  one  could  correct  the  series  solution  for  the  overestimation  caused 
by  the  heat  generation.  For  one-dimensional  heat  transfer,  the  disturbance  in  the 
region  of  the  heated  disk  would  be  composed  of  two  separate  disturbances,  the 
surface  heat  generation  and  the  external  disturbance.  Superposition  holds  because 
the  equations  are  linear.  Therefore,  the  resulting  surface  temperature  changes  axe 
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the  sum  of  the  changes  due  to  the  two  separate  disturbances.  Given  the  total  sur¬ 
face  temperature  changes,  the  series  solution  should  accurately  estimate  the  total 
surface  heat  flux  and  could  be  corrected  to  yield  the  heat  flux  due  to  the  external 
disturbance  alone  by  subtracting  the  value  of  the  surface  heat  generation.  For 
longer  times,  one  can  not  make  this  correction  because  the  heat  transfer  from  the 
localized  surface  heat  generation  becomes  two-dimensioned.  Even  though  super¬ 
position  still  holds,  the  correction  needed  for  the  series  solution  will  be  something 
less  than  the  value  of  the  surface  heat  generation. 

Varying  the  parameter  L/Rcyl  does  not  affect  the  results.  Figure  25  shows 
results  with  the  parameter  L/Rcyl  varying  using  a  value  of  10.0  for  /?,  1/8  for 
Rd,sk / Rcyli  0-1  for  the  tolerance  and  values  of  0.05  and  0.2  for  7.  The  percent 
difference  is  small  and  less  than  the  expected  accuracy  of  the  model.  Figure  65  in 
Appendix  C.l  shows  similar  results  using  a  value  of  100.0  for  f3. 

Varying  the  parameter  RD1SK/RCYL  also  does  not  affect  the  results.  Figure  26 
shows  the  results  with  the  parameter  RD,SK/RCYL  varying  using  a  value  of  10.0 
for  /?,  1/8  for  RDISK  /  RCYL,  0.1  for  the  tolerance  and  values  of  0.05  and  0.2  for  7. 
The  percent  difference  is  again  very  small  and  less  than  the  expected  error  in  the 
model.  Figure  66  in  Appendix  C.l  shows  similar  results  using  a  value  of  100.0  for 
/?.  The  percent  differences  are  larger  but  still  less  than  the  expected  accuracy  of 
the  model. 

4-2  Non- Adiabatic  Cases  with  No  Heat  Generation  or  Preheating 

The  non-adiabatic  cases  with  no  heat  generation  or  preheating  investigate 
the  effect  on  the  series  solution  of  the  two-dimensional  heat  transfer  which  results 
from  allowing  heat  flux  across  the  outer  radial  boundary  of  the  cylinder  in  the 
absence  of  heat  generation  or  preheating.  The  preheating  problem  is  not  used  for 
these  cases,  and  the  parameter  7  is  set  to  zero. 

The  effects  of  allowing  heat  flux  across  the  outer-radial  boundary  of  the 


cylinder  are  bounded  by  the  effects  seen  in  the  limiting  cases.  Each  of  the  three 
limiting  cases  for  heat  flux  across  the  outer- radial  boundary  of  the  cylinder  is 
implemented  using  the  single-material  problem.  The  limiting  case  for  heat  flux 
out  across  the  outer- radial  boundary  of  the  cylinder  uses  Equation  (35)  for  the 
outer-radial  boundary  condition  on  the  cylinder.  The  limiting  case  for  heat  flux  in 
across  the  outer-radial  boundary  of  the  cylinder  uses  Equation  (36)  for  the  outer- 
radial  boundary  condition  on  the  cylinder.  The  case  for  no  heat  flux  across  the 
outer-radial  boundary  of  the  cylinder  is  the  adiabatic  case. 

In  the  absence  of  heat  generation  and  preheating,  the  heat  transfer  in  the 
adiabatic  case  will  be  one-dimensional,  so  the  series  solution  should  be  the  true 
solution  for  times  up  to  the  time  when  the  leading  edge  of  the  external  disturbance 
reaches  the  back  of  the  cylinder.  The  only  approximation  in  the  series  solution  for 
the  adiabatic  case  with  no  heat  generation  or  preheating  is  caused  by  assuming 
the  surface  temperature  to  be  a  piece-wise  linear  function  of  time  (see  Section  2.4). 
Figures  23  and  63  show  that  using  a  value  of  0.0  for  7  in  the  adiabatic  case  does 
indeed  yield  a  percent  difference  very  near  0.0.  Thus,  the  percent  difference  due 
to  any  intermediate  condition  causing  heat  flux  out  across  the  outer-radial  bound¬ 
ary  of  the  cylinder  will  be  bounded  by  zero  and  the  percent  difference  found  in 
the  limiting  case  for  heat  flux  out.  Likewise,  the  percent  difference  due  to  any 
intermediate  condition  causing  heat  flux  in  across  the  outer-radial  boundary  of  the 
cylinder  will  be  bounded  by  zero  and  the  percent  difference  found  in  the  limiting 
case  for  heat  flux  in. 

Two  intermediate  cases  for  heat  flux  across  the  outer-radial  boundary  of  the 
cylinder  are  also  investigated.  Both  use  the  two-material  problem  with  a  geometry 
parameter  RMAX  / RCyl  e(lual  to  2.0.  The  intermediate  case  for  heat  flux  out  models 
teflon  as  the  insulating  material  in  a  test  specimen.  The  property  ratios  k'/k  and 
( Pcp)/(pcp )'  are  -33  and  .49,  respectively,  [7],  [4:p.68S]  and  [6:p.608].  The  boundary 
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condition  on  the  outer- radius  of  the  surrounding  material  for  this  case  is 

6+(r+  =  2.0)  =  0  (121) 

The  intermediate  case  for  heat  flux  in  across  the  boundary  uses  property  ratios  of 
1.0  and  6.0  for  k'Jk  and  (/?Cp)/(pCp)',  respectively.  The  boundary  condition  on  the 
outer-radius  of  the  surrounding  material  for  this  case  is 

d0+\ 

^L„-°  (l22) 

Because  the  non-adiabatic  cases  with  no  heat  generation  or  preheating  do 
not  use  the  preheating  problem,  the  parameter  /?  does  not  really  apply.  Instead, 
values  of  0.1  and  1.0  are  used  for  the  parameter  ,{BiFHONT,  the  Biot  number  at  the 
front  surface.  The  geometry  parameter  L/Rcyl  is  varied  between  0.2  and  5.0,  and 
the  geometry  parameter  RDISK / RCyl  *s  varied  between  the  three  values  of  1/6,  1/8 
and  1/12. 

Figures  27,  28,  29  and  30  show  the  results  for  the  limiting  case  for  heat  flux 
out,  the  intermediate  case  for  heat  flux  out,  the  limiting  case  for  heat  flux  in  and  the 
intermediate  case  for  heat  flux  in,  respectively,  as  the  geometry  parameter  L/RCYl 
varies  using  a  value  of  10.0  for  the  parameter  (3  and  1/8  for  the  geometry  parameter 
Rdisk / RCyl-  As  seen  in  the  results,  the  series  solution  underestimates  the  external 
surface  heat  flux  for  the  cases  with  outward  radial  heat  flux  and  overestimates  the 
external  surface  heat  flux  for  cases  with  inward  radial  heat  flux.  For  all  cases, 
there  is  a  bounding  value  for  the  geometry  parameter  L/Rcyl  such  that  any  lower 
value  for  L/Rcyl  gives  a  minimal  percent  difference.  For  the  limiting  case  for  heat 
flux  out,  limiting  L/Rcyl  to  1.0  keeps  the  percent  difference  less  than  2.0.  For 
the  intermediate  case  for  heat  flux  out,  limiting  L/Rcyl  to  2.0  keeps  the  percent 
difference  less  than  6.0.  For  the  limiting  case  for  heat  flux  in,  limiting  L/Rcyl  to 
0.6  keeps  the  percent  difference  less  than  5.0.  Finally,  for  the  intermediate  case  for 
heat  flux  in,  limiting  L/Rcyl  to  1.0  keeps  the  percent  differnce  less  than  2.0. 
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Figure  27.  Results  for  the  Limiting  Case  for  Heat  Flux  Out  Across  the  Outer- 
Radial  Boundary  with  No  Heat  Generation  as  L/Rcyl  Varies  Using 
front  = 
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Figure  2S.  Results  for  the  Intermediate  Case  for  Heat  Flux  Out  Across  the  Outer- 
Radial  Boundary  with  No  Heat  Generation  as  L/Rcyl  Varies  Using 
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Figure  29.  Results  for  the  Limiting  Case  for  Heat  Flux  In  Across  the  Outer- 
Radial  Boundary  with  No  Heat  Generation  as  L/Rcyl  Varies  Using 
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Figure  30.  Results  for  the  Intermediate  Case  for  Heat  Flux  In  Across  the  Outer- 
Radial  Boundary  with  No  Heat  Generation  as  L/RcyL  Varies  Using 
front  = 


For  values  of  L/Rcyl  larger  than  these  bounding  values,  the  percent  differ¬ 
ences  become  large.  For  the  limiting  case  for  heat  heat  flux  out,  using  a  value 
of  2.0  for  L[Rcyl  and  0.1  for  dBiFRONT  gives  a  steadily  increasing  percent  differ¬ 
ence  which  reaches  30  by  non-dimensional  time  of  0.07.  For  the  intermediate  case 
for  heat  flux  out,  using  a  value  of  4.0  for  L/Rcyl  and  0.1  for  dBiFRONT  gives  a 
steadily  increasing  percent  difference  which  reaches  13  by  non-dimensional  time  of 
0.07.  For  the  limiting  case  for  heat  flux  in,  using  a  value  of  0.8  for  L/Rcyl  and 
0.1  for  dBimosT  fPves  a  rapidly  increasing  percent  difference  which  reaches  85  by 
non-dimensional  time  of  0.07.  Finally,  for  the  intermediate  case  for  heat  flux  in, 
using  a  value  of  2.0  for  L/Rcyl  and  0.1  for  dBiFRONT  gives  a  steadily  increasing 
percent  difference  which  reaches  34  by  non-dimensional  time  of  0.07.  In  all  cases, 
the  error  in  the  series  solution  increases  with  time,  so  one  may  obtain  better  accu¬ 
racy  with  larger  values  for  the  parameter  L/Rcyl  if  one  is  only  interested  in  results 
for  shorter  times.  Although  the  results  from  the  limiting  case  for  heat  flux  in  are 
severe,  one  should  note  that  the  conditions  which  this  limiting  case  model  are  quite 
severe  and  unlikely  to  approximate  any  actual  conditions  in  the  laboratory. 

For  a  short  time  after  the  disturbance,  the  heat  transfer  in  the  region  of 
the  disk  should  be  approximately  one-dimensional.  Therefore,  the  series  solution 
should  theoretically  be  accurate  for  short  times.  Some  of  the  results  clearly  show 
this  delay  in  the  onset  of  error  in  the  series  solution.  As  would  be  expected,  the 
delay  in  terms  of  non-dimensional  time  is  greater  with  smaller  values  of  L/Rcvl 
since  the  tendency  to  produce  radial  gradients  is  less.  As  the  value  of  L/Rcyl 
decreases,  the  resistance  to  heat  transfer  in  the  radial  direction  increases  relative 
to  the  resistance  to  heat  transfer  in  the  axial  direction,  so  the  radial  gradients 
should  be  smaller  when  using  smaller  values  for  L/Rcyl. 

The  trend  for  the  series  solution  to  overestimate  or  underestimate  depend¬ 
ing  on  the  direction  of  the  establishing  radial  gradients  makes  sense  and  can  be 
explained  by  examining  typical  flux  plots  for  the  transient  heat  transfer  in  the 
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cylinder  [4:pp. 135-137].  A  flux  plot  is  a  network  of  isotherms  and  heat  flow  lines. 
The  heat  flow  lines  are  drawn  with  arrows  to  indicate  the  direction  of  the  heat 
flow  and  temperature  gradient.  Heat  flow  lines  must  always  be  perpendicular  to 
isotherms.  The  area  between  adjacent  heat  flow  lines  is  termed  a  lane.  Heat  can 
be  visualized  as  flowing  in  these  lanes.  In  a  well  drawn  flux  plot,  the  heat  flux 
between  two  adjacent  isotherms  in  a  lane  can  be  estimated  by 

AT 

qA,a  *  *  Ad  (123) 

where  q  is  the  heat  flux,  k  is  the  thermal  conductivity,  At  is  the  approximate  area 
of  the  lane  perpendicular  to  the  heat  flow  line,  AT  is  the  temperature  difference 
between  the  isotherms  and  Ad  is  the  approximate  distance  between  the  isotherms. 
Heat  energy  must  be  conserved,  so  any  heat  which  does  not  continue  to  flow  in  the 
lane  will  cause  the  temperature  to  rise.  Figure  31  shows  a  typical  flux  plot  for  heat 
transfer  in  an  initially  isothermal,  semi-infinite  solid.  In  contrast,  Figure  32  shows 
typical  flux  plots  for  heat  transfer  in  the  intermediate  cases  for  heat  flux  out  and 
heat  flux  in  across  the  outer-radial  boundary. 

The  series  solution  uses  the  time  history  of  changes  in  the  surface  temperature 
to  estimate  the  surface  heat  flux.  The  series  solution  is  accurate  only  if  the  changes 
in  surface  temperature  that  do  occur  are  equal  to  those  that  would  occur  in  an 
initially  isothermal,  semi-infinite  solid  with  the  same  time  history  of  surface  heat 
flux.  For  the  cases  that  establish  outward  radial  heat  flux,  the  changes  in  surface 
temperature  for  any  given  time  history  of  surface  heat  flux  are  less  than  the  changes 
in  surface  temperature  that  would  occur  in  an  initally  isothermal,  semi-infinite 
solid  with  the  same  time  history  of  surface  heat  flux.  Given  the  smaller  changes 
in  surface  temperature,  the  series  solution  underestimates  the  surface  heat  flux. 

The  trend  for  surface  temperature  changes  to  be  less  in  runs  that  establish 
outward  radial  gradients  can  be  visualized  by  examining  the  typical  flux  plots. 
Figures  31  and  32  show  flux  plots  for  heat  transfer  in  a  semi-infinite  solid  and 
in  the  cylinder  with  outward  radial  gradients,  respectively.  As  can  be  seen  by 
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Figure  31.  Typical  Flux  Plot  for  Heat  Transfer  in  an  Initially  Isothermal,  Semi- 
Infinite  Solid 
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a.  Heat  Flux  Out  Across  the  Outer-Radial  Boundary 


b.  Heat  Flux  In  Across  the  Outer-Radial  Boundary 


Figure  32.  Typical  Flux  Plots  for  Heat  Transfer  in  the  Intermediate  Cases  for 
Heat  Flux  Out  and  In  Across  the  Outer- Radial  Boundary  of  the  Cylin¬ 
der 


the  increasing  area  of  the  lanes  in  its  flux  plot,  the  heat  flux  within  the  solid 
with  outward  radial  gradients  will  be  less  than  the  heat  flux  within  the  semi¬ 
infinite  solid  when  the  same  surface  heat  flux  is  present  in  both  in  accordance 
with  Equation  (123).  Then,  the  temperature  changes  which  result  from  the  same 
time  history  of  surface  heat  flux  will  be  less  since  smaller  temperature  changes 
are  required  to  establish  the  heat  flux  through  the  cylinder.  With  inward  radial 
gradients,  the  area  of  the  lanes  decreases,  so  the  temperature  changes  will  be 
greater,  and  the  series  solution  overestimates. 

Figures  67,  68,  69  and  70  in  Appendix  C.2  show  the  results  for  the  limiting 
case  for  heat  flux  out,  the  intermediate  case  for  heat  flux  out,  the  limiting  case  for 
heat  flux  in  and  the  intermediate  case  for  heat  flux  in  as  the  parameter  L/Rcyl 
varies  using  a  value  of  1.0  for  dBiFROKT  instead.  The  results  using  a  larger  value 
for  dBiFRONT  show  the  very  same  trends  but  with  slightly  less  percent  differences. 
It  must  be  noted,  however,  that  the  dimensional  error  in  the  series  solution  is  still 
greater  for  the  larger  value  of  dBiFRONT  since  the  percent  difference  is  normalized 
with  the  value  of  dBiFRONT  in  the  denominator.  The  smaller  percent  difference 
when  using  the  larger  value  for  dBiFRONT  may  be  misleading  also  because  of  the 
way  the  results  are  displayed.  It  is  likely  that  the  fractional  amount  of  the  true 
surface  heat  flux  by  which  the  series  solution  overestimates  or  underestimates  stays 
fairly  constant  when  the  value  for  dBiFRONT  varies.  The  way  to  show  this  would  be 
to  evaluate  a  new  percent  difference  using  the  finite-difference  estimate  for  surface 
heat  flux  as  the  normalizing  factor.  The  non-dimensional  temperatures  at  the 
front  surface  in  the  disturbance  problem  will  tend  to  increase  when  using  larger 
values  for  dBiFRONT ,  so  non-dimensional  surface  heat  flux  values  decrease.  If,  in 
fact,  the  series  solution  errs  by  a  constant  fractional  amount  of  the  finite-difference 
estimate  when  the  value  for  dBiFRONT  changes,  then  the  percent  difference  is  less 
for  larger  values  of  dBiFRONT  only  because  the  finite- difference  estimates  for  surface 
heat  fluxes  are  smaller. 
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Figure  33.  Results  for  the  Limiting  Case  for  Heat  Flux  Out  Across  the  Outer 
Radial  Boundary  with  No  Heat  Generation  as  RDtSK/RCYL  Varies  Us 
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Figure  35.  Results  for  the  Limiting  Case  for  Heat  Flux  In  Across  the  Outer- 
Radial  Boundary  with  No  Heat  Generation  as  RDISK/ RCYl  Varies  Us¬ 
ing  dBi FRONT  ~ 
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Figure  36.  Results  for  the  Intermediate  Case  for  Heat  Flux  In  Across  the  Outer- 
Radial  Boundary  with  No  Heat  Generation  as  RDISK / Varies  Us¬ 
ing  jBi  front  —  •  ^ 


The  geometry  parameter  RDISK/RCYL  has  little  impact  on  the  results.  Fig¬ 
ures  33,  34,  35  and  36  show  the  results  as  the  parameter  RDISk/^cyl  varies  using 
a  value  of  0.1  for  dBiFRONT-  Varying  RDISk/RCyl  *n  the  absence  of  heat  generation 
merely  changes  the  area  over  which  the  surface  temperature  is  averaged  for  use  in 
the  series  solution.  In  all  cases  except  the  intermediate  case  for  heat  flux  out,  a 
larger  value  for  RDISk/Rcyl  produces  slightly  more  error  in  the  series  solution.  This 
trend  makes  sense  since  the  heat  transfer  should  be  more  one-dimensional  closer  to 
the  centerline  of  the  cylinder  and  since  the  radial  derivative  is  zero  directly  at  the 
centerline.  The  impact  of  varying  the  geometry  parameter  RDISk/Rcyl  increases 
tor  larger  values  of  L/Rcyl  which  also  makes  sense  since  the  radial  derivatives  are 
greater  with  larger  values  of  L/Rcyl.  However,  it  should  be  noted  that  the  change 
in  the  percent  difference  as  the  parameter  L/Rcyl  varies  in  most  cases  is  less  than 
the  expected  accuracy  in  the  model,  so  the  results  are  inconclusive.  The  trends 
could  be  trends  in  the  error  equation  of  the  finite-difference  model  rather  than 
actual  trends  in  the  results.  Figures  71,  72,  73  and  74  in  Appendix  C.2  show  the 
results  as  Rdisk/Rcyl  varies  using  a  value  of  1.0  for  dBiFRONT.  These  results  are 
very  similar. 

4-3  Non- Adiabatic  Casts  with  Heat  Generation  and  Preheating 

The  non-adiabatic  cases  with  heat  generation  and  preheating  investigate  the 
effect  on  the  series  solution  when  both  causes  of  two-dimensional  heat  transfer 
are  present.  These  runs  are  a  subset  of  the  runs  in  the  previous  section  to  which 
preheating  and  subsequent  heat  generation  are  added. 

Both  the  preheating  and  the  disturbance  problems  are  used  for  all  runs.  It 
should  be  noted  that  the  limiting  cases  for  heat  flux  out  and  in  are  identical  in  the 
preheating  problem  since  the  outer-radial  boundary  condition  is  the  same  in  both 
as  is  seen  by  Equations  27  and  28  in  Section  2.2.  The  limiting  case  for  heat  flux  in 
models  the  condition  where  the  thermal  conductivity  in  the  surrounding  material 
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is  infinite  and  the  temperature  of  the  surrounding  material  is  the  temperature  of 
the  fluid  at  the  front  surface.  In  the  preheating  problem,  the  temperature  of  the 
fluid  at  the  front  surface  is  the  initial  temperature. 

The  tolerance  used  in  the  convergence  criterion  for  the  preheating  problem  is 
0.1  for  all  runs.  All  runs  use  a  value  of  1/8  for  the  geometry  parameter  RDISK/RCYL  • 
Each  of  the  four  cases  is  rim  with  a  few  values  for  the  geometry  parameter  L/Rcyl. 
Values  of  0.0,  0.05  and  0.2  for  the  parameter  7  axe  run  with  each  of  the  values  for 
L/Rcyl.  All  runs  are  duplicated  using  values  of  10.0  and  100.0  for  the  parameter 
0. 

Figures  37,  38,  39  and  40  show  the  results  for  the  limiting  case  for  heat  flux 
out,  the  intermediate  case  for  heat  flux  out,  the  limiting  case  for  heat  flux  in  and 
the  intermediate  case  for  heat  flux  in,  respectively,  using  a  value  of  10.0  for  /?. 
Figures  75,  76,  77  and  78  in  Appendix  C.3  show  the  results  for  the  limiting  case 
for  heat  flux  out,  the  intermediate  case  for  heat  flux  out,  the  limiting  case  for  heat 
flux  in  and  the  intermediate  case  for  heat  flux  in,  respectively,  using  a  value  of  100.0 
for  /?.  As  in  the  adiabatic  cases,  the  preheating  and  subsequent  heat  generation 
has  no  impact  on  the  results. 
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Figure  37.  Results  for  the  Limiting  Case  for  Heat  Flux  Out  Across  the  Outer 
Radial  Boundary  with  Heat  Generation  and  Preheating  Using  = 
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Figure  3S. 


Results  for  the  Intermediate  Case  for  Heat  Flux  Out  Across  the  Outer- 
Radial  Boundary  with  Heat  Generation  and  Preheating  Using  = 
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Figure  39.  Results  for  the  Limiting  Case  for  Heat  Flux  In  Across  the  Outer- Radial 
Boundary  with  Heat  Generation  and  Preheating  Using  f3  =  10.0 
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Figure  40.  Results  for  the  Intermediate  Case  for  Heat  Flux  In  Across  the  Outer- 
Radial  Boundary  with  Heat  Generation  and  Preheating  Using  j3  = 


V.  Conclusions  and  Recommendations 


This  thesis  numerically  investigates  the  accuracy  of  the  one-dimensional  series 
solution  in  determining  the  external  surface  heat  flux  actually  present  at  the  film 
in  thin-film  heat  transfer  gages.  The  heat  transfer  problem  is  simplified  to  two 
dimensions.  From  an  overall  point  of  view,  the  two  possible  causes  for  error  in 
the  series  solution  are  the  electrical  heat  generation  in  the  thin-film  sensor  and  the 
radial  heat  transfer  induced  by  non-adiabatic  conditions  at  the  boundary  of  the 
gage. 

The  results  show  that  the  electrical  heat  generation  should  not  cause  errors 
as  long  as  the  instrumentation  is  turned  on  and  the  transients  from  the  electrical 
heat  generation  are  allowed  to  settle  down  prior  to  the  test.  It  should  not  take 
long  for  the  transients  to  settle  down  sufficiently  to  keep  the  error  small.  In  tests 
using  an  adiabatic  condition  at  the  boundary  of  the  gage,  a  non-dimensional  time 
of  0.118,  which  is  approximately  twice  as  long  as  the  time  it  takes  the  leading 
edge  of  a  thermal  disturbance  to  travel  through  the  gage,  was  sufficient.  The 
transients  should  settle  down  even  quicker  with  non-adiabatic  conditions  since  the 
temperature  changes  will  be  less. 

The  radial  heat  transfer  induced  from  non-adiabatic  conditions  at  the  bound¬ 
ary  of  the  gage  can  cause  very  significant  errors  in  the  results  for  long  times.  Be¬ 
cause  of  limitations  on  the  model,  this  investigation  only  looks  at  non-dimensional 
times  greater  than  0.01.  The  two  most  significant  non-dimensional  parameters 
which  influence  the  error  caused  by  the  non-adiabatic  conditions  are  the  ratio  of 
the  thermal  diffusivities  in  the  insulating  and  cylinder  materials  and  the  geome¬ 
try  parameter  L/Rcyl.  Outward  radial  gradients  caused  by  using  an  insulating 
material  with  a  smaller  thermal  diffusivity  than  that  of  the  cylinder  material  will 
cause  the  measured  heat  flux  to  be  less  than  the  actual,  external  surface  heat  flux. 
Inward  radial  gradients  caused  by  using  an  insulating  material  with  a  greater  ther- 
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mal  diffusivity  will  cause  the  measured  heat  flux  to  be  greater  than  the  actual.  In 
either  case,  a  larger  value  of  the  geometry  parameter  L/Rcyl  increases  the  error. 

Limiting  cases  which  bound  the  effects  caused  by  non-adiabatic  conditions 
were  investigated.  In  the  limiting  case  for  outward  radial  heat  flux,  keeping  the 
geometry  parameter  L/Rcyl  less  than  1.0  limited  the  error  in  the  measurement 
for  heat  flux  to  less  than  2  percent  of  the  theoretically  maximum,  external  surface 
heat  flux  for  non-dimensional  times  up  to  0.07.  In  the  limiting  case  for  inward 
radial  heat  flux,  keeping  the  geometry  parameter  L/Rcyl  less  than  0.6  limited  the 
error  to  5  percent.  In  an  intermediate  case  for  outward  radial  heat  flux,  which 
used  a  value  of  1/6  for  the  ratio  of  the  thermal  diffusivities  of  the  surrounding 
material  and  the  cylinder  material,  keeping  the  geometry  parameter  L/Rcyl  less 
than  2.0  limited  the  error  to  6  percent.  Finally,  in  an  intermediate  case  for  inward 
radial  heat  flux,  which  used  a  ratio  of  6.0  for  the  thermal  diffusivities,  keeping  the 
geometry  parameter  L/Rcyl  less  than  1.0  limited  the  error  to  2  percent. 

Using  values  of  L/Rcyl  larger  than  these  can  produce  large  errors  in  the 
results.  In  the  limiting  case  for  outward  radial  heat  flux,  using  a  value  of  2.0  for 
the  geometry  parameter  L/Rcyl  and  0.1  for  the  Biot  number  at  the  front  surface 
produced  a  steadily  increasing  error  which  reached  30  percent  by  non-dimensional 
time  of  0.07.  In  the  limiting  case  for  inward  radial  heat  flux,  using  a  value  of 
0.8  for  the  geometry  parameter  L/Rcyl  and  0.1  for  the  Biot  number  at  the  front 
surface  produced  an  error  which  reached  85  percent.  In  the  intermediate  case  for 
outward  radial  heat  flux,  using  a  value  of  4.0  for  L/Rcyl  and  0.1  for  the  Biot 
number  produced  an  error  which  reached  13  percent.  Finally,  in  the  intermediate 
case  for  inward  radial  heat  flux,  using  a  value  of  2.0  for  L/Rcyl  and  0.1  for  the  Biot 
number  produced  an  error  which  reached  34  percent.  Since  the  error  increases  with 
time,  one  can  attain  better  accuracy  with  larger  values  of  the  geometry  parameter 
L/Rcyl  and  greater  differences  in  the  thermal  diffusivities  if  one  uses  shorter  test 
times. 
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The  error  caused  by  a  non-adiabatic  condition  at  the  boundary  of  the  gage 
can  be  significant.  One  can  try  to  limit  this  error  in  three  ways.  First,  one  can  use 
as  good  an  insulating  material  as  possible  in  mounting  the  gage.  Although  this 
technique  promises  to  improve  the  accuracy  in  the  series  solution,  it  may  cause  the 
heat  flux  at  the  gage  to  be  very  different  from  the  heat  flux  that  would  occur  in  the 
absence  of  the  gage  because  the  temperature  changes  due  to  the  electrical  heating 
will  be  greater.  Secondly,  one  could  try  to  match  properties  between  the  cylinder 
material  and  the  insulating  material  to  limit  the  temperature  difference  across  the 
boundary  of  the  cylinder.  This  technique  should  improve  the  accuracy  of  the  series 
solution  while  minimizing  the  deleterious  effect  of  the  electrical  heating  in  the  thin 
film.  Thirdly,  one  can  use  a  gage  with  a  smaller  geometry  ratio  L/Rcyl. 

The  results  from  this  investigation  must  be  weighed  in  view  of  the  accuracy 
of  the  model  as  estimated  from  the  transient  check  cases  for  one-dimensional  heat 
transfer  in  the  social  direction.  As  a  conservative  estimate,  the  model  should  be 
accurate  to  within  a  percent  difference  of  4.0  for  non-dimensional  times  greater 
than  0.01  where  the  percent  difference  is  the  percent  of  the  theoretically  maximum, 
external  surface  heat  flux  by  which  the  series  solution  estimate  differs  from  the 
actual  heat  flux  in  the  finite-difference  model,  i.e.  (g+£/8  ~  d  )  x  100.00  . 
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Appendix  A.  Nodal  Equations 


A. I  Single- Material  Problem  with  Fully- Insulated,  Outer-Radial  Boundary 
Condition 


Preheating  Problem: 

Note  1:  Nodal  Equations  (a)-(i)  apply  to  the  corresponding  positions  in 
Figure  4. 

Note  2:  For  nodes  on  the  front  surface  of  the  cylinder,  Equations  (a)-(c),  the 
last  term  in  the  nodal  finite-difference  equation  is  included  only  if  the  node  lies 
within  the  heated  disk. 
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Disturbance  Problem: 

Note  1:  Nodal  Equations  (a)-(i)  apply  to  the  corresponding  positions  in 
Figure  4. 

Note  2:  For  nodes  on  the  front  surface  of  the  cylinder,  Equations  (a)-(c), 
the  last  term  in  the  nodal  finite-difference  equation,  i.e.  the  term  including  the 
parameter  7,  is  included  only  if  the  node  lies  within  the  heated  disk. 
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A. 2  Two-Material  Problem  with  Pally-Insulated,  Outer-Radial  Boundary 
Condition 


Preheating  Problem: 

Note  1:  Nodal  Equations  (a)-(o)  apply  to  the  corresponding  positions  in 
Figure  5. 

Note  2:  For  nodes  on  the  front  surface  of  the  cylinder,  Equations  (a)-(c),  the 
last  term  in  the  nodal  finite-difference  equation  is  included  only  if  the  node  lies 
within  the  heated  disk. 
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Disturbance  Problem; 


Note  1:  Nodal  Equations  (a)-(o)  apply  to  the  corresponding  positions  in 
Figure  5. 

Note  2:  For  nodes  on  the  front  surface  of  the  cylinder,  Equations  (a)-(c), 
the  last  term  in  the  nodal  finite-difference  equation,  i.e.  the  term  including  the 
parameter  7,  is  included  only  if  the  node  lies  within  the  heated  disk. 
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Appendix  B.  Graphs  of  Results  from  the  Check  Cases 


Steady-State  Check  Cases  for  Heat  Transfer  «'n  the  Axial  Direction 
B.1.1  Preheating  Problems 


Prohoatlng  Problem: 


132 


Figure  41.  Other  Results  for  Steady  State  Heat  Transfer  in  the  Axial  Direction 
Using  the  Preheating  Problem,  Figure  (a) 


Preheating  Problem 


Figure  42.  Other  Results  for  Steady  State  Heat  Transfer  in  the  Axial  Direction 
Using  the  Preheating  Problem,  Figure  (b) 


B.l.S  Disturbance  Problems 
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Disturbance  Problem 


Figure  43.  Other  Results  for  Steady  State  Heat  Transfer  in  the  Axial  Direction 
Using  the  Disturbance  Problem,  Figure  (a) 


FRONT 


Results  for  Steady  State  Heat  Transfer  in  the  Axial  Direction 
the  Disturbance  Problem,  Figure  (b) 


Disturbance  Problem:  6+  vs. 
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Figure  45.  Other  Results  for  Steady  State  Heat  Transfer  in  the  Axial  Direction 
Using  the  Disturbance  Problem,  Figure  (c) 
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'Ults  for  Steady  State  Heat  Transfer  in  the  Axial  Direction 
Disturbance  Problem,  Figure  (d) 


Disturbance  Problem 


Results  for  Steady  State  Heat  Transfer  in  the  Axial  Direction 
the  Disturbance  Problem,  Figure  fe) 


B.i  Steady- State  Cheek  Cases  for  Heat  Transfer  in  the  Radial  Direction 

Case  numbers  refer  to  the  cases  given  in  Tables  2  and  3.  See  Figure  13  in 
Section  3.2  for  the  results  from  Case  2  using  the  single-material  problem. 

B.i.l  Single- Material  Problems 
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Figure  48.  Case  1  for  Steady  State  Heat  Transfer  in  the  Radial  Direction  Using 
the  Single-Material  Problem 
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Single-Material  Problem: 


Figure  49.  Case  3  for  Steady  State  Heat  Transfer  in  the  Radial  Direction  Using 
the  Single-Material  Problem 
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Case  4: 


B.l.t  Two- Material  Problems 
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Case  2: 

Boundary  Conditions 
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Figure  53.  Case  3  for  Steady  State  Heat  Transfer  in  the  Radial  Direction  Using 
the  Two-Material  Problem 


Figure  54.  Case  4  for  Steady  State  Heat  Transfer  in  the  Radial  Direction  Using 
the  Two- Material  Problem 
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B.S  Transient  Check  Cases 

B.S.l  The  Semi-Infinite  Solid  with  Constant  Surface  Heat  Flux 
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Finite-Difference  Model 


t+ 


Figure  55.  Finite-Difference  Model  Values  for  6+  vs.  t+  for  the  Semi-Infinite  Solid 
with  Constant  Surface  Heat  Flux 


Series  Estimate  for  qf+ 

(Using  the  Finite-Difference  Estimates  for  9f+)  vs.  t+ 


t+ 


Figure  57.  Series  Solution  Estimates  for  qf  (Using  the  Finite-Difference  Model 
Values  for  Of)  vs.  t+  for  the  Semi-Infinite  Solid  with  Constant  Surface 
Heat  Flux 


B.S.t  The  Plane  Wall  with  Convection 
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Figure  5S.  Finite- Difference  Model  Values  for  6+  vs.  t+  for  the  Plane  Wall  with 
Convection 
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Figure  61.  Series  Solution  Estimates  for  q +  (Using  the  Finite-Difference  Model 
Values  for  )  vs.  t+  for  the  Plane  Wall  with  Convection 


<+ 


Figure  62.  Comparing  the  Finite- Difference  Estimates.  Series  Solution  Estimates 
(Using  the  Finite- Difference  Model  Values  for  $+)  and  Analytical  Val¬ 
ues  for  q+  vs.  t+  for  the  Plane  Wall  with  Convection 
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Figure  64.  Results  with  Heat  Generation  and  Preheating  as  the  Tolerance  Varies 
Using  0  =  100.0 
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Results  for  the  Intermediate  Case  for  Heat  Flux  Out  Across  the  Outer- 
Radial  t  undary  with  No  Heat  Generation  as  L/Rcyl  Varies  Using 
d^i  F  vr  1.0 
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Figure  69.  Results  for  the  Limiting  Case  for  Heat  Flux  In  Across  the  Outer- 
Radial  Boundary  with  No  Heat  Generation  as  L/Rcyl  Varies  Using 

d^i  FRONT 
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Q+.o)  *  100.00  vs.  £+ 


Figure  72.  Results  for  the  Intermediate  Case  for  Heat  Flux  Out  Across  the  Outer- 
Radial  Boundary  with  No  Heat  Generation  as  RDtsK/ RCYL  Varies  Us¬ 
ing  dBiFROST  =  1-0 


Figure  73.  Results  for  the  Limiting  Case  for  Heat  Flux  In  Across  the  Outer- 
Radial  Boundary  with  No  Heat  Generation  as  RDISl</ RCyl  Varies  Us 
ing  dRi front  ~  UO 


9t„.)  x  100  00 


t+ 


Figure  75.  Results  for  the  Limiting  Case  for  Heat  Flux  Out  Across  the  Outer- 
Radial  Boundary  with  Heat  Generation  and  Preheating  Using  ( 3  = 
100.0 
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Figure  77.  Results  for  the  Limiting  Case  for  Heat  Flux  In  Across  the  Outer- Radial 
Boundary  with  Heat  Generation  and  Preheating  Using  [S  =  100.0 


176 


)  x  100.00 


Figure  78. 


Results  far  the  Intermediate  Case  for  Heat  Flux  In  Across  the  Outer- 
Radial  Boundary  with  Heat  Generation  and  Preheating  Using  jS  = 


100.0 
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Appendix  D.  Example  Programs 


D.l  The  Finite-Difference  Model  for  the  Single- Material  Problem  with 


C 

C 

C 

C 

C 

C 

C 
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c 


c 

c 


c 
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c 

c 

c 

c 
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Fully -Insulated,  Outer-Radial  Boundary  Condition 
PROGRAM  GAGE2H 

************************************************** 

*  2LT  JOSEPH  A.  BONAFEDE,  GA-88M 

*  FALL  1987 

* 

*  ADVISOR:  DR.  JAMES  E.  HITCHCOCK 
**************************************************, 


GLOSSARY  OF  MAIN  VARIABLES: 

THAOLD  -  ARRAY  STORING  NON-DIMENSIONAL  TEMPERATURE  VALUES  FOR 
ALL  GRID  POINTS.  USE  THAOLD  TO  STORE  INITIAL  VALUES 
AND  TO  STORE  VALUES  DURING  PREVIOS  TIME  STEP. 

USED  INTERNALLY 

THANEW  -  ARRAY  TO  STORE  NEW  NON-DIMENSIONAL  TEMPERATURE  VALUES 
CALCULATED  DURING  A  NEW  TIME  STEP. 

USED  INTERNALLY 

IHAX  -  DETERMINES  NUMBER  OF  NODES  IN  THE  RADIAL  DIRECTION.  NUMBER 
OF  NODES  IN  THE  RADIAL  DIRECTION  EQUALS  (IMAX  ♦  1). 

INPUT  PARAMETER 

KMAX  -  DETERMINES  THE  NUMBER  OF  NODES  IN  THE  Z  DIRECTION  (FRONT  TO 
BACK  OF  GAGE) .  NUMBER  OF  NODES  IN  THE  Z  DIRECTION  EQUALS 
(KMAX  +1). 

INPUT  PARAMETER 

IGEN  -  ARRAY  WHICH  STORES  INFORMATION  ABOUT  TOP  SURFACE  NODES. 

STORES  VALUE  OF  1.0  FOR  A  NODE  IF  THE  NODE  IS  IN  THE  HEAT 
GENERATING  REGION  AND  A  VALUE  OF  0.0  FOR  A  NODE  IF  IT  IS  NOT. 
USED  INTERNALLY 

IGENMX  -  LARGEST  RADIAL  NODE  INCLUDED  IN  THE  HEAT  GENERATING  REGION. 
INPUT  PARAMETER 

DELT  -  TIME  STEP  (NON-DIMENSIONAL  TIME). 

INPUT  PARAMETER 

MAXT  -  MAXIMUM  NUMBER  OF  TIME  STEPS  ALLOWED. 

INPUT  PARAMETER 

TDIST  -  TIME  STEP  AT  WHICH  TO  CHANGE  TO  THE  EXTERNAL  DISTURBANCE 
PROBLEM  (AS  LONG  AS  TDIST  .LE.  MAXT). 

INPUT  PARAMETER 

TIME1  -  TIME  STEP  AT  WHICH  TO  CHANGE  TO  THE  EXTERNAL  DISTURBANCE 
PROBLEM  IN  THE  PROGRAM. 
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C  TIME1-0  *>  EXTERNAL  DISTURBANCE  PROBLEM  ALWAYS 

C  2<»TIME1<«MAXT  «>  SWITCH  FROM  NON-DISTURBANCE  PROBLEM 

C  TO  DISTURBANCE  PROBLEM  AT  N-TIME1 

C  TIME1*MAXT+1  *>  NON-DISTURBANCE  PROBLEM  ALWAYS 

C  USED  INTERNALLY 

C  COUNT  -  NUMBER  OF  TIME  STEPS  INTO  THE  DISTURBANCE  PROBLEM  (I.E. 

C  THE  NUMBER  OF  TIME  STEPS  ALREADY  RUN  FOR  THE  DISTURBANCE 

C  PROBLEM) . 

C  USED  INTERNALLY 

C  MAXTHA  -  LARGEST  ABSOLUTE  VALUE  FOR  NEW  NON-DIMENSIONAL  TEMPERATURE 
C  IN  THE  TIME  STEP  (I.E.  LARGEST  VALUE  STORED  IT  THANEW 

C  ARRAY) . 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


USED  INTERNALLY 

CHNG  -  ABSOLUTE  CHANGE  IN  NON-DIMENSIONAL  TEMPERATURE  AT  EACH 
NODE  DURING  THE  LAST  TIME  STEP. 

USED  INTERNALLY 

VISCHNG  -  CHANGE  IN  NON-DIMENSIONAL  TEMPERATURE  AT  EACH  NODE  DURING 
THE  LAST  TIME  STEP  AS  A  FRACTION  OF  THE  LARGEST  VALUE  FOR 
NON-DIMENSIONAL  TEMPERATURE  IN  THAT  TIME  STEP  (I.E. 

VISCHNG  *  CHNG/MAXTHA  ) . 

USED  INTERNALLY 

MAXCH  -  MAXIMUM  VALUE  FOR  VISCHNG  DURING  THE  TIME  STEP  (AS  LONG  AS 

THE  VALUE  FOR  MAXTHA  FOR  THE  TIME  STEP  IS  NOT  NEAR  ZERO.  IF 
THE  VALUE  FOR  MAXTHA  IS  NEAR  ZERO.  USE  ABSOLUTE  CHANGE  RATHER 
THAN  VISUAL  CHANGE  AS  THE  COVERGENCE  CRITERIA.) 

USED  INTERNALLY 

TOL  -  VALUE  TO  USE  IN  DETERMINING  CONVERGENCE  (STEADY  STATE  PROBLEMS 
ONLY).  IF  RATE  OF  VISUAL  CHANGE  IS  LESS  THAN  TOL,  THEN  ASSUME 
THAT  THE  PROGRAM  HAS  CONVERGED  SUFFICIENTLY. 

(I.E.  MAXCH/DELT  <  TOL  *>  CONVERGENCE  BECAUSE  RATE  OF  CHANGE 

IS  SUFFICIENTLY  SMALL  ) 

INPUT  PARAMETER 

LRRAT  -  THE  GEOMETRY  RATIO  (LENGTH  OF  GAGE) /(RADIUS  OF  GAGE). 

INPUT  PARAMETER 

GBIOT  -  THE  TOP  SURFACE  BIOT  NUMBER  BEFORE  THE  DISTURBANCE. 

INPUT  PARAMETER 

GBIOTB  -  THE  BOTTOM  SURFACE  BIOT  NUMBER  BEFORE  THE  DISTRUBANCE. 

INPUT  PARAMETER 

DBIOT  -  THE  TOP  SURFACE  BIOT  NUMBER  AFTER  THE  DISTRUBANCE. 

INPUT  PARAMETER 

DBIOTB  -  THE  BOTTOM  SURFACE  BIOT  NUMBER  AFTER  THE  DISTURBANCE. 

INPUT  PARAMETER 

BIOT  -  TOP  SURFACE  BIOT  NUMBER  AT  THE  PARTICULAR  TIME  STEP. 

USED  INTERNALLY 
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C  BIOTB  -  BOTTOM  SURFACE  BIOT  HUMBER  AT  THE  PARTICULAR  TIME  STEP. 

C  USED  INTERNALLY 

C  GENRAT  -  THE  “GENERATION  RATIO"  (QGEN/DBIOT)/(TFFINAL  -  TINIT) . 

C  INPUT  PARAMETER 

C  BFRAT  -  THE  RATIO  (TFBACK  -  TINIT) /(TFFRONT  -  TINIT)  WHERE  TFBACK 

C  AND  TFFRONT  ARE  THE  FLUID  TEMPERATURES  DURING  THE 

C  DISTURBANCE  PROBLEM. 

c  input  parameter 

C  DGBIRAT  -  THE  RATIO  (DBIOT)/(GBIOT) . 

C  USED  INTERNALLY 

C  DELR  -  SPACIAL  STEP  IN  THE  NON-DIMENSIONAL  RADIAL  DIRECTION. 

C  USED  INTERNALLY 

C  DELZ  -  SPACIAL  STEP  IN  THE  NON-DIMENSIONAL  Z  DIRECTION. 

C  USED  INTERNALLY 

Cl-  INDEX  VARIABLE  FOR  NON-DIMENSIONAL  RADIAL  DIRECTION. 

C  USED  INTERNALLY 

C  K  -  INDEX  VARIABLE  FOR  NON-DIMENSIONAL  Z  DIRECTION. 

C  USED  INTERNALLY 

C  N  -  INDEX  VARIABLE  FOR  TIME  STEPS. 

C  USED  INTERNALLY 

C  RUNNUM  -  THE  RUN  NUMBER  (USED  FOR  BOOK  KEEPING  PURPOSES)  . 

INPUT  VARIABLE 

:  LAKR1  ,LAMZ1  ,LAMZ2  ,P1  -  COMMON  PRODUCT  TERMS  IN  THE  FINITE 

ELEMENT  EQUATIONS. 

USED  INTERNALLY 

DECLARE  VARIABLES: 

IMPLICIT  CHARACTER(A-Z) 

REAL  THAOLDCO : 20 ,0 : 20) .THANEW (0:20,0:20) 

REAL  IGEN(0:20) 

REAL  DELT 

INTEGER  MAXT.TDIST.TIMEl. COUNT 
REAL  MAXTHA , CHNG , VISCHNG , MAXCH , TOL 
INTEGER  IMAX , KMAX , IGENMX 
REAL  LRRAT 

REAL  BIOT, BIOTB, GBIOT,GBIOTB,DBIOT,DBIOTB 
REAL  DGBIRAT, GENRAT, BFRAT 
REAL  DELR, DELZ, LAMR1.LAMZ1.LAMZ2, PI 
REAL  A1,A2,A3,A,B,C1,C,D 
INTEGER  RUNNUM 
INTEGER  I.K.N 
INTEGER  X.Y 
C 

C . — . 
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C  OPEM  INPUT/OUTPUT  FILES. 

C 

OPEN  (UNIT-1  .FILE- ’G2INP  * .  STATUS  - ' OLD » ) 

OPEN  (UNIT-2  .FILE- * G20UT » .  STATUS-  *  NEW ' ) 

OPEN (UNIT-3  .FILE- » INI X  * , STATUS- ’ OLD » ) 

OPEN  (UNIT- 10 ,  FILE-  *  DISKT »  .STATUS-  *  NEW  * ) 

REWIND(UNIT-l) 

REWIND(UNIT«3) 

C 

C . . 

READ  INPUT  PARAMETERS. 

READ(1,*) 

READ(1,*)  DELT 
READ(1,*)  MAXT 
READ(1,*)  TDIST 
READ(1,*)  TOL 
C 

READ(1,*)  IMAX 
READ(1,*)  KMAX 
READ(1,*)  IGENMX 
READ(1,*)  LRRAT 

C 

READ(1,*)  GBIOT 
READ(1,*)  GBIOTB 
READ(1,*)  DBIOT 
READ(1.*)  DBIOTB 
READ(1,*)  GENRAT 
READ(1,*)  BFRAT 
READ(1,*)  DGBIRAT 
READ(1,*)  RUNNUM 
C 

C . 

C  INITIALIZE  PARAMETER  NEEDED  IN  PRINTING  OUT  HEADER  (  USING  HEADER 
C  INSTEAD  OF  JUST  ECHOING  THE  INPUT  )  . 

C 

IF  (TDIST  .LE.  1  )  THEN 
TIME1  -  0 

ELSEIF  (TDIST  .GT.  MAXT)  THEN 
TIME1  «  MAXT  ♦  1 
ELSE 

TIME1  -  TDIST 
ENDIF 
C 
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C  PRINT  HEADER. 

C 

WRITE(2,*)  'RUN  NUMBER  \RUNNUM 
HRITEC2.*) 'INPUT  VALUES:' 

HRITE(2,*) 

WRITEC2,*) ’**  GEOMETRY  OF  THE  GAGE  AND  FINITE  ELEMENT  MESH  **’ 
WRITEC2,*) 

WRITE(2,*) ’  GRID  POINTS  IN  R-DIRECTION  :  0  -  ’.IMAX 
WRITEC2,*)’  GRID  POINTS  IN  Z-DIRECTION  :  0  -  * ,KMAX 
IF  (IGENMX  .LT.  0)  THEN 

WRITE(2,*) ’  R-DIRECTION  GRID  POINTS  INCLUDED’ 

HRITEC2,*)’  IN  THE  HEAT  GENERATING  DISK  :  NONE’ 

ELSEIF  (IGENMX  .EQ.  0)  THEN 

WRITE(2,*)’  R-DIRECTION  GRID  POINT  INCLUDED’ 

WRITE(2,*)’  IN  THE  HEAT  GENERATING  DISK  :  0’ 

ELSE 

WRITE(2,*)’  R-DIRECTION  GRID  POINTS  INCLUDED* 

WRITE (2,*) ’  IN  THE  HEAT  GENERATING  DISK  :  0  -  ’.IGENMX 
ENDIF 
WRITEC2,*) 

WRITE(2,*) ’  RATIO  OF  (LENGTH  OF  GAGE) /(RADIUS  OF  GAGE)  =  * , LRRAT 
WRITE(2,*) 

WRITE(2,*) 

WRITE(2 ,  *)  ’ **TIME  STEP  AND  LIMITS  INCLUDING  EXTERNAL  PARAMETERS**’ 
WRITE(2,*) 

WRITE(2.*)’  DELTA  T  *  ’.DELT 
WRITE(2,*) 

WRITE(2,*) ’  MAXIMUM  NUMBER  OF  TIME  STEPS  *  ’ ,MAXT 
WRITE(2,*) ’  TOLERANCE  FOR  CONVERGENCE  «  ’ ,TOL 
WRITE(2,*) 

IF  (TIME1  .EQ.  2)  THEN 

WRITE(2,*) ’  TIME  STEP  ASSIGNED  TO  THE  "START  UP"’ 

WRITE(2,*)  ’  PROBLEM  (  NO  EXTERNAL  DISTURBANCE  )  :  1» 

WRITE(2,*) ’  BIOT  NUMBER  AT  THE  TOP  SURFACE  »  ’  .GBIOT 

WRITE(2,*) ’  BIOT  NUMBER  AT  THE  BOTTOM  SURFACE  =  ’ .GBIOTB 

ENDIF 

IF  (TIME1  .GT.  2)  THEN 

WRITE (2,*)’  TIME  STEPS  ASSIGNED  TO  THE  "START  UP"’ 

WRITE(2,*) ’  PROBLEM  (  NO  EXTERNAL  DISTURBANCE  )  :  1  -  ’. 

t  TIME1-1 

WRITE(2,*)  ’  BIOT  NUMBER  AT  THE  TOP  SURFACE  *  ’.GBIOT 

WRITE(2,*) ’  BIOT  NUMBER  AT  THE  BOTTOM  SURFACE  *  ’  .GBIOTB 

ENDIF 
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IF  ((TIMEl  .EQ.  0) .AND. (MAXT  .GT.  D)  THEM 
WRITE (2,*) 

WRITE(2,*)’  TIME  STEPS  ASSIGNED  TO  THE  PROBLEM' 

WRITEC2.*)'  WITH  AN  EXTERNAL  DISTURBANCE  :  1  -  MAXT 
WRITE(2,*)’  BIOT  NUMBER  AT  THE  TOP  SURFACE  «  ’ .DBIOT 

WRITE(2,*) '  BIOT  NUMBER  AT  THE  BOTTOM  SURFACE  ■  ’.DBIOTB 

WRITEC2.*)’  RATIO  (QGEN/DBIOT)/(TFFIMAL  -  TI)  »  ’.GENRAT 

WRITEC2,*)'  RATIO  (TFBACK  -  TI)/(TFFROMT  -  TI)  -  ’.BFRAT 

WRITEC2,*)’  RATIO  (DBIOT)/(GBIOT)  ■  \DGBIRAT 

ENDIF 

IF  ( (TIMEl  .EQ.  0). AND. (MAXT  .EQ.  1))  THEN 
WRITEC2,*) 

WRITEC2,*)’  TIME  STEP  ASSIGNED  TO  THE  PROBLEM’ 

WRITEC2,*)’  WITH  AN  EXTERNAL  DISTURBANCE  :  1* 

WRITE(2,*)’  BIOT  NUMBER  AT  THE  TOP  SURFACE  =  DBIOT 

WRITE(2,*) ’  BIOT  NUMBER  AT  THE  BOTTOM  SURFACE  =  ’.DBIOTB 

WRITEC2,*)’  RATIO  (QGEN/DBIOT) / (TFFINAL  -  TI)  =  ’ .GENRAT 

WRITE(2,*) ’  RATIO  (TFBACK  -  TI)/ (TFFRONT  -  TI)  =  ’.BFRAT 

WRITE(2,*) ’  RATIO  (DBIOT)/(GBIOT)  =  ’ .DGBIRAT 

ENDIF 

IF  ((TIMEl  .GT.  0). AND. (TIMEl  .EQ.  MAXT))  THEN 
WRITE (2,*) 

WRITE(2,«0  ’  TIME  STEP  ASSIGNED  TO  THE  PROBLEM’ 

WRITE (2,*) ’  WITH  AN  EXTERNAL  DISTURBANCE  :  ’.TIMEl 

WRITE(2,*)  ’  BIOT  NUMBER  AT  THE  TOP  SURFACE  *  ’.DBIOT 

WRITE(2,*) ’  BIOT  NUMBER  AT  THE  BOTTOM  SURFACE  *  ’.DBIOTB 

WRITE(2,*) ’  RATIO  (QGEN/DBIOT) /(TFFINAL  -  TI)  =  ’.GENRAT 

WRITE(2,*) ’  RATIO  (TFBACK  -  TI)/ (TFFRONT  -  TI)  =  ’.BFRAT 

WRITE(2 ,*) ’  RATIO  (DBIOT)/ (GBIOT)  *  ’.DGBIRAT 

ENDIF 

IF  ((TIMEl  .GT.  0). AND. (TIMEl  .LT.  MAXT))  THEN 
WRITE(2,*) 

WRITE(2,*) ’  TIME  STEPS  ASSIGNED  TO  THE  PROBLEM’ 

WRITE(2,*) ’  WITH  AN  EXTERNAL  DISTURBANCE  :  ’, 
t  TIMEl,’  -  ’.MAXT 

WRITE(2,*) ’  BIOT  NUMBER  AT  THE  TOP  SURFACE  =  ’.DBIOT 

WRITE(2  ,*)  ’  BIOT  NUMBER  AT  THE  BOTTOM  SURFACE  =  ’.DBIOTB 

WRITE(2 ,*) ’  RATIO  (QGEN/DBIOT) /(TFFINAL  -  TI)  =  ’.GENRAT 

WRITE(2,*) ’  RATIO  (TFBACK  -  TI)/ (TFFRONT  -  TI)  =  ’.BFRAT 

WRITE(2,*) ’  RATIO  (DBIOT)/(GBIOT)  *  ’.DGBIRAT 

ENDIF 
WRITE (2.*) 

WRITE(2,*) 

WRITEC2.*) ’RESULTS:’ 
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WRITEC2,*) 


C  INITIALIZE  ARRAYS  AND  OTHER  PARAMETERS. 

C 

READC3,*) 

DO  100  K  -  O.KMAX 
DO  110  I  «  O.IMAX 
THANEW(I.K)  *  0. 

READ(3,5)  THAOLD(I.K) 

IF  (TIME1  .EQ.  0)  THAOLD(I.K)  «  THAOLD(I ,K)*GENRAT*DGBIRAT 
110  CONTINUE 

100  CONTINUE 
C 

DO  120  I  «  O.IMAX 

IF  (I  .LE.  IGENMX)  THEN 
IGEN(I)  *  1. 

ELSE 

IGEN(I)  =  0. 

END  IF 
120  CONTINUE 
C 

DELR  «  1. /FLOAT (IMAX) 

DELZ  «  1 . /FLOAT(KMAX) 

C 

LAMR1  ■  DELT/DELR**2 
LAMZ1  ■  DELT/DELZ**2 
LAMZ2  *  DELT/DELZ 
C 

PI  =  (FLOAT(IMAX)  -  .5)/(FL0AT(IMAX)  -  .25) 

C 

. . 

c  ************************************************************** 

C  . 

C  MAIN  LOOP. 

C 


IF  CTIME1  .EQ.  0)  THEN 
BIOT  *  DBIOT 
BIOTB  *  DBIOTB 
ELSE 

BIOT  *  GBIOT 
BIOTB  *  GBIOTB 
END  IF 
COUNT  »  0 
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MAXCH  ■  99.99 
C 

. . 

DO  500  N  «  l.MAXT 
IF  ((MAXCH/DELT)  .GT.  TOL)  THEM 
C 
C 
C 
C 
C 
C 
C 


1010 
1000 

C 
C 

C  SOLVE  FOR  NEXT  TIME  STEP. 

C 

C 

C  FRONT  OF  GAGE 

C 
C 
C 

C  (A) 

C 

A  *  4.*(LRRAT**2)*LAMR1*  THAOLD(l.O) 

B  «  2. *LAMZ1*  THAOLD(O.l) 

Cl  *  1.  -  4.*(LRRAT**2)*LAMR1  -  2.*LAMZ1  -  2.*LAMZ2*BI0T 
C  *  Cl*  THAOLD(O.O) 

IF  (N  .LT.  TIMED  THEN 

D  ■  2.*LAMZ2*BI0T*  IGEN(O) 

ELSE 

D  *  2.*LAMZ2*BI0T*(1 .  ♦  IGEN(0)*GENRAT) 

ENDIF 

THANEW(O.O)  -A+B+C+D 
C 

C  (B) 


REDEFINE  BIOT  NUMBERS  AND  UPDATE  THAOLD  WHEN  STARTING 
THE  DISTURBANCE  PROBLEM. 

IF  (N  .EQ.  TIMED  THEN 
BIOT  *  DBIOT 
BIOTB  *  DBIOTB 
DO  iOOO  I  *  O.IMAX 

DO  1010  K  »  O.KMAX 

THAOLD(I.K)  «  THAOLD(I , K) *GENRAT*DGBIRAT 
CONTINUE 
CONTINUE 
ENDIF 
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DO  600  I  -  l.IMAX-1 

*1-1.  -  l./(2.*FL0AT(I)) 

*2  -  1.  ♦  i./(2.*FL0AT(I)) 

A3  ■  (LRRAT**2) *LAMR1 

A  -  A3*(  Al*  THAOLD(I-l.O)  ♦  A2*  THA0LD(I+1 .0)  ) 

B  -  2.*LAMZ1*  THAOLDd , 1) 

Cl  -  1.  -  2 . * (LRRAT**2) *LAMR1  -  2.*LAMZ1  -  2 . *LAMZ2*BI0T 
C  »  Cl*  THAOLD(I ,0) 

IF  (N  .LT.  TIMED  THEN 

D  «  2.*LAMZ2*BI0T*  IGEN(I) 

ELSE 

D  «  2.*LAMZ2*BI0T*(1 .  *  IGEN(I)*GENRAT) 

ENDIF 

THANEW(I ,0)  -A+B+C+D 
600  CONTINUE 

C 

C  (C) 

C 

A  *  2. *(LRRAT**2)*LAMR1*P1*  THA0LD(IMAX-1 ,0) 

B  «  2. *LAMZ1*  THAOLD(IMAX, 1) 

Cl  *  1.  -  2.*(LRRAT**2)*LAMR1*P1  -  2.*LAMZ1  -  2. *LAMZ2*BI0T 
C  «  Cl*  THA0LD(IMAX,0) 

IF  (N  .LT.  TIMED  THEN 

D  *  2.*LAMZ2*BI0T*  IGEN(IMAX) 

ELSE 

D  «  2 . *LAMZ2*BI0T* ( 1 .  *  I GEN (I MAX ) *GENRAT ) 

ENDIF 

THANEW(IMAX.O)  -A+B+C+D 
C 
C 

C  INTERIOR  OF  GAGE 

C 

C 

DO  700  K  *  l.KMAX-l 
C 

C  (D) 

C 

A  «  4.*(’  1RAT**2)*LAMR1*  THA0LD(1 ,K) 

B  «  LAMZ1*(  THAOLD(O.K-l)  ♦  THA0LD(0,K+1)  ) 

Cl  -  1.  *  4 . *(LRRAT**2) *LAMR1  -  2.*LAMZ1 
C  »  Cl*  THA0LD(0,lO 
THANEN(O.K)  •  A  *  B  ♦  C 
C 

C  (E) 
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C 


DO  750  I  •  l.IMAX-l 

A1  «  1.  -  l./(2.*FL0AT(I)> 

A2  •  1.  ♦  l./(2.*FL0AT(I)) 

A3  •  (LRRAT**2) *LAMR1 

A  -  A3*<  Al*  THA0UKI-1  ,K)  ♦  A2*  THA0LD(I*1,K)  ) 
B  -  LAMZ1*(  THA0LD(I,K-1)  ♦  THA0LD(I,K+1)  ) 

Cl  «  1.  -  2 . * (LRRAT**2) *LAMR1  -  2.*LANZ1 
C  «  Cl*  THAOLDCI.K) 

THAHEW(I,K)  «  A  ♦  B  ♦  C 
7S0  CONTINUE 

C 

C  (F) 

C 

A  *  2  •  * (LRRAT **2) *LAMR 1 *P 1 *  THAOLD (IMAX- 1 , K) 

B  *  LAMZ1*(  THAOLD ( IKAX ,  K- 1 )  +  THAOLD(IMAX ,K+1)  ) 

Cl  *  1.  -  2.*(LRRAT**2)*LAMR1*P1  -  2.*LAMZ1 
C  *  Cl*  THAOLD (IHAX.K) 

THANEW(IMAX.K)  *  A  ♦  B  ♦  C 
C 

700  CONTINUE 


C  BACK  OF  GAGE 

C 
C 
C 

C  (G) 

C 

A  »  4.*(LRRAT**2)*LAHR1*  THA0LD(1 ,KMAX) 

B  «  2.*LAKZ1*  THAOLD (O.KMAX-1) 

Cl  «  1.  -  4.*(LRRAT**2)*LAMR1  -  2.*LAMZ1  -  2. *LAMZ2*BI0TB 
C  *  Cl*  THAOLD(O.KHAX) 

IF  (N  .LT.  TIMED  THEN 
D  *  0- 
ELSE 

0  *  2 . *LAMZ2*BI0TB*  BFRAT 
ENDIF 

THANEW(O.KMAX)  *A+B+C*D 
C 

C  CH) 

C 

DO  800  I  *  l.IMAX-l 

Al  *  1.  -  l./(2.*FL0AT(I>) 
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A2  «  1.  ♦  l./(2.*FL0AT(I)) 

A3  »  (LRRAT**2)*LAMR1 

A  -  A3*(  Al*  THAOLD (1-1. KMAX)  ♦  A2*  THAOLO (!♦ 1 , KMAX )  ) 

B  ■  2.*LAMZ1*  THAOLD (I .KMAX-l) 

Cl  «  1.  -  2 . * (LRRAT**2) ♦LAMR1  -  2.*LAMZ1  -  2.*LAMZ2*BI0TB 
C  -  Cl*  THAOLD(I.KMAX) 

IF  (M  .LT.  TIMED  THEN 


D  »  0. 

ELSE 

D  ■  2 . *LAMZ2*BI0TB*  BFRAT 
ENDIF 

THANEW(I.KMAX)  »A  +  B  +  C*D 
800  CONTINUE 

C 

C  (I) 

C 

A  =  2.*(LRRAT**2)*LAMR1*P1*  THA0LD(IMAX-1 .KMAX) 

B  =  2.*LAMZ1*  THAOLD(IMAX, KMAX-l) 

Cl  *  1.  -  2.*(LRRAT**2)*LAMR1*P1  -  2.*LAMZ1  -  2.*LAMZ2*BI0TB 
C  =  Cl*  THAOLD(IMAX.KMAX) 

IF  (N  .LT.  TIMED  THEN 
D  ■  0. 

ELSE 

D  «  2 . *LAMZ2*BI0TB*  BFRAT 
ENDIF 

THANEW(IMAX.KMAX)  -A+B+C+D 
C 

C  . 

C  WRITE  NON-DIMENSIONAL  TEMPERATURE  VALUES  AT  THE  HEATED  DISK 

C  TO  FILE  "DISKT"  (DURING  DISTURBANCE  PROBLEM  ONLY). 

C  **  MODIFIED  TO  WRITE  TO  "DISKT"  ONLY  AT  INITIAL  TIME  STEP 

C  OF  DISTURBANCE  PROBLEM  AND  AT  TIME  STEPS  T*  *  .001-. 070 

C  WITH  DELT-.001  AFTER  THE  DISTURBANCE  PROBLEM  STARTS. 

C  BE  CAREFUL  TO  HAVE  ENOUGH  TIME  STEPS  TO  HAVE  ATLEAST 

C  .070  IN  NON-DIMENSIONAL  TIME  FOR  THE  DISTURBANCE  (I.E. 

C  NEED  TO  HAVE  COUNT  ATLEAST  EQUAL  TO  70  WHEN  THE  PROGRAM 

C  FINISHES). 

C 

IF  (N  .GE.  TIMED  THEN 

IF  (  (N.EQ. TIMED  -OR.  ((TIME1 .EQ .0) . AND. (N.EQ . 1) )  )  THEN 
COUNT  «  0 

WRITE (10, 25)  COUNT 
DO  1100  I  ■  O.IGENMX 

WRITE(10,30)  I.THAOLD(I.O) 
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1100  CONTINUE 

ENDIF 
C 

X  «  N  -  MAX(TIMEl.l)  ♦  1 
Y  »  NINT ( . 001/DELT) 

IF  (  MOD(X.Y)  .EQ.  0  )  THEN 
COUNT  *  COUNT  ♦  1 
WRITE (10,25)  COUNT 
DO  1125  I  «  O.IGENMX 

WRITE(10,30)  I.THANEW(I.O) 

1125  CONTINUE 

ENDIF 
ENDIF 

C 

C  . 

c  ROLL  DOWN  THAOLD(1:20,1:20)  AND  FIND  MAXIMUM  VISUAL 

C  PERCENT  CHANGE.  (BE  CAREFUL  NOT  TO  DIVIDE  BY  ZERO.) 

C 

MAXTHA  «  0. 

DO  900  I  *  O.IMAX 
DO  910  K  *  O.KMAX 

IF  (  ABS(THANEW(I,K))  .GT.  MAXTHA)  THEN 
MAXTHA  «  ABS (THANEW ( I , K ) ) 

ENDIF 

910  CONTINUE 

900  CONTINUE 

C 

MAXCH  »  0. 

DO  950  I  *  O.IMAX 
DO  960  K  *  O.KMAX 

CHNG  «  ABS(  THANEW(I.K)  -  THAOLD(I.K)  ) 

IF  (  MAXTHA  .GT.  l.E-15)  THEN 
VISCHNG  ■  CHNG/MAXTHA 
IF  (VISCHNG  .GT.  MAXCH)  MAXCH  *  VISCHNG 
ELSE 

IF  (CHNG  .GT.  MAXCH)  MAXCH  «  CHNG 
ENDIF 

HAOLD(I.K)  *  THANEW ( I, K) 

960  CONTINUE 

950  CONTINUE 

C 

C . 

ELSE 

C  . 
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C 

CONVERGES.  JUMP  OUT  OF  LOOP. 

WRITE(2,*) 

WRITE(2,*) ’CONVERGES  AFTER  ITERATIONS.' 

WRITE  (2,*) 'MAXIMUM  CHANGE  ON  THE  LAST  ITERATION  WAS  ’.MAXCH.’.’ 
GOTO  2000 
C 

. . 

ENDIF 

SOO  CONTINUE 

. . 

C 

C  IF  YOU  EXECUTE  THESE  STATEMENTS,  THEN  YOU  COMPLETED  THE 

C  MAXIMUM  NUMBER  OF  ITERATIONS  FOR  THE  LOOP  WITHOUT  CONVERGING. 

C  PRINT  MESSAGE  IF  EXPECTING  SOME  CONVERGENCE  (I.E.  TOL  >  0). 

C 

IF  (TOL  .GT.  0.)  THEN 
WRITE(2,*) 

WRITE(2,*) 'QUITS  WITHOUT  CONVERGING  AFTER  \N-1,’  ITERATIONS.’ 
WRITEC2.*) 'MAXIMUM  CHANGE  ON  THE  LAST  ITERATION  WAS  ’  ,MAXCH, ’ .  ’ 
ENDIF 
C 

. . 

c  ******************************************************************* 

. . 

2000  CONTINUE 
C 

. . - . 

C  PRINT  CHECK  ON  VALUE  OF  COUNT. 

C 

IF  (COUNT  .LT.  70)  THEN 
WRITE(2,*) 

WRITE(2,*) 'WARNING:  COUNT  LESS  THAN  70.  COUNT  =  ’.COUNT 
WRITE(2,*) 

ELSE 

WRITE(2,*) 

WRITE(2,*) 'VALUE  OF  COUNT  IS  ’.COUNT 
WRITE(2,*) 

WRITE(2,*) 

ENDIF 

C 

. . 

C  FORMAT  STATEMENTS. 
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F0RMAT(20X,E17. 10) 
25  FORKAT(IS) 

30  F0RMAT(I5,5X,E17 . 10) 

C 

END 
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D.t  The  Finite-Difference  Model  for  the  Two-Material  Problem  with 
Fully -Insulated,  Outer-Radial  Boundary  Condition 
PROGRAM  GAGE4M 

C  *************************************************** 

C  *  2LT  JOSEPH  A.  BONAFEDE,  GA-88M  * 

C  *  FALL  1987  * 

C  *  ♦ 

C  *  ADVISOR:  DR.  JAMES  E.  HITCHCOCK  * 

C  *************************************************** 

c 

c  GLOSSARY  OF  MAIN  VARIABLES: 

C  THAOLD  -  ARRAY  STORING  NON-DIMENSIONAL  TEMPERATURE  VALUES  FOR 
C  ALL  GRID  POINTS.  USE  THAOLD  TO  STORE  INITIAL  VALUES 

C  AND  TO  STORE  VALUES  DURING  PREVIOUS  TIME  STEP. 

C  USED  INTERNALLY 

C  THANEW  -  ARRAY  TO  STORE  NEW  NON-DIMENSIONAL  TEMPERATURE  VALUES 
C  CALCULATED  DURING  A  NEW  TIME  STEP. 

C  USED  INTERNALLY 

C  IMAX  -  DETERMINES  NUMBER  OF  NODES  IN  THE  RADIAL  DIRECTION.  NUMBER 
C  OF  NODES  IN  THE  RADIAL  DIRECTION  EQUALS  (IMAX  ♦  1) . 

C  INPUT  PARAMETER 

C  KMAX  -  DETERMINES  THE  NUMBER  OF  NODES  IN  THE  Z  DIRECTION  (FRONT  TO 
C  BACK  OF  GAGE) .  NUMBER  OF  NODES  IN  THE  Z  DIRECTION  EQUALS 

C  (KMAX  ♦  1). 

C  INPUT  PARAMETER 

C  IGAGE  -  DETERMINES  THE  NUMBER  OF  NODES  IN  THE  RADIAL  DIRECTION 
C  WHICH  ARE  INCLUDED  IN  THE  GAGE.  THE  RADIAL  DIRECTION 

C  NODES  WHICH  MAKE  UP  THE  GAGE  ARE  NODES  WITH  1  =  0-  IGAGE. 

C  ALSO  DETERMINES  THE  NON-DIMENSIONAL  STEP  SIZE  IN  THE 

C  RADIAL  DIRECTION.  THE  NON-DIMENSIONAL  STEP  SIZE  IN  THE 

C  RADIAL  DIRECTION  EQUALS  1. /IGAGE  (I.E.  DELR  *  1. /IGAGE). 

INPUT  PARAMETER 

IGEN  -  ARRAY  WHICH  STORES  INFORMATION  ABOUT  TOP  SURFACE  NODES. 

STORES  VALUE  OF  1.0  FOR  A  NODE  IF  THE  NODE  IS  IN  THE  HEAT 
GENERATING  REGION  AND  A  VALUE  OF  0.0  FOR  A  NODE  IF  IT  IS  NOT. 
USED  INTERNALLY 

IGENMX  -  LARGEST  RADIAL  NODE  INCLUDED  IN  THE  HEAT  GENERATING  REGION. 
C  INPUT  PARAMETER 

C  DELT  -  TIME  STEP  (NON-DIMENSIONAL  TIME) . 

C  INPUT  PARAMETER 

C  MAXT  -  MAXIMUM  NUMBER  OF  TIME  STEPS  ALLOWED. 

C  INPUT  PARAMETER 

C  TDIST  -  TIME  STEP  AT  WHICH  TO  CHANGE  TO  THE  EXTERNAL  DISTURBANCE 

C  PROBLEM  (AS  LONG  AS  TDIST  .LE.  MAXT). 
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INPUT  PARAMETER 

TIME1  -  TIME  STEP  AT  WHICH  TO  CHANGE  TO  THE  EXTERNAL  DISTURBANCE 
PROBLEM  IN  THE  PROGRAM. 

TIME1-0  ■>  EXTERNAL  DISTURBANCE  PROBLEM  ALWAYS 
2<*TIME1<*MAXT  «>  SWITCH  FROM  NON-DISTURBANCE  PROBLEM 
TO  DISTURBANCE  PROBLEM  AT  H-TIME1 
TIME1»MAXT*1  «>  NON-DISTURBANCE  PROBLEM  ALWAYS 
USED  INTERNALLY 

COUNT  -  NUMBER  OF  TIME  STEPS  INTO  THE  DISTURBANCE  PROBLEM  (I.E. 

THE  NUMBER  OF  TIME  STEPS  ALREADY  RUN  FOR  THE  DISTURBANCE 
PROBLEM) . 

USED  INTERNALLY 

MAXTHA  -  LARGEST  ABSOLUTE  VALUE  FOR  NEW  NON-DIMENSIONAL  TEMPERATURE 
IN  THE  TIME  STEP  (I.E.  LARGEST  VALUE  STORED  IN  THANEW 
ARRAY) . 

USED  INTERNALLY 

CHNG  -  ABSOLUTE  CHANGE  IN  NON-DIMENSIONAL  TEMPERATURE  AT  EACH 
NODE  DURING  THE  LAST  TIME  STEP. 

USED  INTERNALLY 

VISCHNG  -  CHANGE  IN  NON-DIMENSIONAL  TEMPERATURE  AT  EACH  NODE  DURING 
THE  LAST  TIME  STEP  AS  A  FRACTION  OF  THE  LARGEST  VALUE  FOR 
NON-DIMENSIONAL  TEMPERATURE  IN  THAT  TIME  STEP  (I.E. 

VISCHNG  «  CHNG/MAXTHA  ). 

USED  INTERNALLY 

MAXCH  -  MAXIMUM  VALUE  FOR  VISCHNG  DURING  THE  TIME  STEP  (AS  LONG  AS 
THE  VALUE  FOR  MAXTHA  FOR  THE  TIME  STEP  IS  NOT  NEAR  ZERO.  IF 
THE  VALUE  FOR  MAXTHA  IS  NEAR  ZERO,  USE  ABSOLUTE  CHANGE  RATHER 
THAN  VISUAL  CHANGE  AS  THE  CONVERGENCE  CRITERIA.) 

USED  INTERNALLY 

TOL  -  VALUE  TO  USE  IN  DETERMINING  CONVERGENCE  (STEADY  STATE  PROBLEMS 
ONLY).  IF  RATE  OF  VISUAL  CHANGE  IS  LESS  THAN  TOL,  THEN  ASSUME 
THAT  THE  PROGRAM  HAS  CONVERGED  SUFFICIENTLY. 

(I.E.  MAXCH/DELT  <  TOL  «>  CONVERGENCE  BECAUSE  RATE  OF  CHANGE 

IS  SUFFICIENTLY  SMALL  ) 

INPUT  PARAMETER 

LRRAT  -  THE  GEOMETRY  RATIO  (LENGTH  OF  GAGE) /(RADIUS  OF  GAGE). 

INPUT  PARAMETER 

GBIOT  -  THE  TOP  SURFACE  BIOT  NUMBER  BEFORE  THE  DISTURBANCE. 

INPUT  PARAMETER 

GBIOTB  -  THE  BOTTOM  SURFACE  BIOT  NUMBER  BEFORE  THE  DISTRUBANCE. 

INPUT  PARAMETER 

DBIOT  -  THE  TOP  SURFACE  BIOT  NUMBER  AFTER  THE  DISTRUBANCE. 

INPUT  PARAMETER 

DBIOTB  -  THE  BOTTOM  SURFACE  BIOT  NUMBER  AFTER  THE  DISTURBANCE. 
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C  IMPUT  PARAMETER 

C  BIOT  -  TOP  SURFACE  BIOT  HUMBER  AT  THE  PARTICULAR  TIME  STEP. 

C  USED  INTERNALLY 

C  BIOTB  -  BOTTOM  SURFACE  BIOT  HUMBER  AT  THE  PARTICULAR  TIME  STEP. 

C  USED  INTERNALLY 

C  GENRAT  -  THE  "GENERATION  RATIO"  (QGEN/DBIOT)/(TFFINAL  -  TINIT) . 

C  INPUT  PARAMETER 

C  BFRAT  -  THE  RATIO  (TFBACK  -  TINIT)/ (TFFRONT  -  TINIT)  WHERE  TFBACK 
C  AND  TFFRONT  ARE  THE  FLUID  TEMPERATURES  DURING  THE 

C  DISTURBANCE  PROBLEM. 

C  INPUT  PARAMETER 

C  DGBIRAT  -  THE  RATIO  (DBIOT)/(GBIOT) . 

C  USED  INTERNALLY 

C  RCGB  -  PROPERTY  RATIO  (RHO*CP)GAGE  /  (RHO*CP)MATERIAL  B. 

C  INPUT  PARAMETER 

C  KBG  -  PROPERTY  RATIO  (K)MATERIAL  B  /  (K)GAGE. 

C  INPUT  PARAMETER 

C  RCGGB  -  PROPERTY  RATIO  (RHO*CP)GAGE/(RHO*CP)WEIGHTED  AVG  GAGE  AND  B. 
C  USED  INTERNALLY 

C  KGBG  -  PROPERTY  RATIO  (K)WEIGHTED  AVG  GAGE  AND  B  /  (K)GAGE. 

C  USED  INTERNALLY 

C  DELR  -  SPACIAL  STEP  IN  THE  NON-DIMENSIONAL  RADIAL  DIRECTION. 

C  USED  INTERNALLY 

C  DELZ  -  SPACIAL  STEP  IN  THE  NON-DIMENSIONAL  Z  DIRECTION. 

C  USED  INTERNALLY 

C  DER  -  VALUE  OF  THE  PARTIAL  DERIVATIVE  IN  THE  RADIAL  DIRECTION 
C  OF  THE  NON-DIMENSIONAL  TEMPERATURE  VARIABLE  AT  THE 

C  OUTER  BOUNDARY  OF  THE  PROBLEM.  USE  A  ONE-SIDED, 

C  BACKWARD  FINITE  DIFFERENCE. 

C  USED  INTERNALLY 

C  MAXDER  -  THE  LARGEST  VALUE  FOR  DER  (IN  ABSOLUTE  VALUE)  DURING 
C  THE  WHOLE  RUN. 

C  USED  INTERNALLY 

Cl-  INDEX  VARIABLE  FOR  NON-DIMENSIONAL  RADIAL  DIRECTION. 

C  USED  INTERNALLY 

C  K  -  INDEX  VARIABLE  FOR  NON-DIMENSIONAL  Z  DIRECTION. 

C  USED  INTERNALLY 

C  N  -  INDEX  VARIABLE  FOR  TIME  STEPS. 

C  USED  INTERNALLY 

C  RUNNUM  -  THE  RUN  NUMBER  (USED  FOR  BOOK  KEEPING  PURPOSES) . 

C  INPUT  VARIABLE 

C  LAMR1 ,LAMZ1 ,LAMZ2,P1 ,P2,P3  -  COMMON  PRODUCT  TERMS  IN  THE  FINITE 
C  ELEMENT  EQUATIONS. 

C  USED  INTERNALLY 
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C  DECLARE  VARIABLES: 

IMPLICIT  CHARACTER(A-Z) 

REAL  THAOLD (0 :  40 , 0 : 40) ,THANEW(0:40,0:40) 
REAL  IGEN(0:20) 

REAL  DELT 

INTEGER  NAXT.TDIST.TIMEl, COUNT 
REAL  MAXTHA , CHNG , VISCHNG , MAXCH , TOL 
INTEGER  IMAX , KM AX , IGENMX , IG AGE 
REAL  LRRAT 

REAL  BIOT,BIOTB,GBIOT,GBIOTB,DBIOT,DBIOTB 

REAL  DGBIRAT , GENRAT , BFRAT 

REAL  RCGB , KBG , RCGGB , KGBG 

REAL  DELR , DELZ , LAMR1 , LAMZ 1 , LAMZ2 , P 1 , P2 , P3 

REAL  DER.MAXDER 

REAL  A1,A2,A3,A,B1,B,C1,C2,C3,C,D 

REAL  NUMER.DENOM 

INTEGER  RUNNUM 

INTEGER  I,K,N 

INTEGER  X.Y 

C 

C . 

C  OPEN  INPUT/OUTPUT  FILES. 

C 

OPEN (UNIT* 1 . FILE* » G4INP » , STATUS* » OLD ’ ) 

OPEN (UNIT=2 . FILE* » G40UT * , STATUS* » NEW ' ) 

OPEN (UNIT*3 , FILE* ' INIT ’ , STATUS* ’ OLD » ) 

OPEN (UNIT*10 , FILE* ’ DISKT 1 , STATUS* ’ NEW  * ) 
REWIND (UNIT*1) 

REV.'IND(UNIT*3) 

C 

C . . 

C  READ  INPUT  PARAMETERS. 

C 

READCl,*) 

READ(1 ,*)  DELT 
READ(1,*)  MAXT 
READCl,*)  TDIST 
READCl,*)  TOL 
C 

READCl,*)  IMAX 
READCl,*)  KMAX 
READCl,*)  IGAGE 
READCl,*)  IGENMX 
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READCl ,*)  LRRAT 
C 

READCl,*)  GBIOT 
READCl,*)  GBIOTB 
READCl,*)  DBIOT 
READCl ,*)  DBIOTB 
READCl ,*)  GEN RAT 
READCl.*)  BFRAT 
READCl,*)  DGBIRAT 
C 

READCl,*)  RCGB 
READCl,*)  KBG 
C 

READCl,*)  RUNNUM 
C 

C . — . 

C  INITIALIZE  PARAMETER  NEEDED  IN  PRINTING  OUT  HEADER  C  USING  HEADER 
C  INSTEAD  OF  JUST  ECHOING  THE  INPUT  ). 

C 

IF  CTDIST  .LE.  1  )  THEN 
TIME1  *  0 

ELSEIF  CTDIST  .GT.  MAXT)  THEN 
TIME1  *  MAXT  +  1 
ELSE 

TIME1  *  TDIST 
ENDIF 

C 

C . 

C  PRINT  HEADER. 

C 

WRITEC2,*) ’RUN  NUMBER  «  ’ , RUNNUM 
WRITEC2,*) 'INPUT  VALUES:’ 

VRITEC2,*) 

HRITEC2,*)’**  GEOMETRY  OF  THE  GAGE  AND  FINITE  ELEMENT  MESH  **’ 
HRITEC2,*) 

WRITEC2,*)’  GRID  POINTS  IN  R-DIRECTION  :  0  -  ’.IMAX 
WRITEC2,*)’  GRID  POINTS  IN  Z-DIRECTION  :  0  -  ’ ,KMAX 
WRITEC2,*)’  R-DIRECTION  GRID  POINTS  INCLUDED’ 

WRITEC2,*)’  IN  THE  GAGE  :  0  -  ’ , IGAGE 

IF  CIGENMX  .LT.  0)  THEN 

HRITEC2,*)’  R-DIRECTION  GRID  POINTS  INCLUDED’ 

HRITEC2,*)’  IN  THE  HEAT  GENERATING  DISK  :  NONE’ 

ELSEIF  CIGENMX  .Eq.  0)  THEN 

WRITEC2,*)'  R-DIRECTION  GRID  POINT  INCLUDED’ 
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IN  THE  HEAT  GENERATING  DISK 


R-DIRECTION  GRID  POINTS  INCLUDED' 
IN  THE  HEAT  GENERATING  DISK 


0  -  '.IGENMX 


' , LRRAT 


WRITEC2,*)' 

ELSE 

WRITE(2,*)’ 

WRITE(2,*) * 

ENDIF 
WRITE(2,*) 

WRITE(2,*)’  RATIO  OF  (LENGTH  OF  GAGE) / (RADIUS  OF  GAGE) 

WRITE(2,*) 

WRITE(2,*) 

WRITE(2 ,  *)  ’ **TIME  STEP  AND  LIMITS  INCLUDING  EXTERNAL  PARAMETERS**’ 
WRITE(2,*) 

WRITE(2,*)  ’  DELTA  T  «  \DELT 
WRITE(2,*) 

HRITE(2,*) ’ 

WRITE(2,*)  ’ 

WRITE(2,*) 

IF  (TIME1  .Eq. 

WRITE(2,*) ’ 

WRITE(2,*) ’ 

WRITE(2,*)’ 

WRITE(2,*) ’ 

WRITE (2,*) 

WRITE(2,*) ’ 

WRITEC2,*)’ 

ENDIF 

IF  (TIME1  .GT.  2)  THEN 

WRITE(2,*) ’  TIME  STEPS  ASSIGNED  TO  THE  "START  UP"’ 

WRITE(2,*) ’  PROBLEM  (  NO  EXTERNAL  DISTURBANCE  )  :  1  - 

t  TIME1-1 

WRITE(2.*)’  BIOT  NUMBER  AT  THE  TOP  SURFACE  =  ’.GBIOT 

WRITE(2,*) ’  BIOT  NUMBER  AT  THE  BOTTOM  SURFACE  =  ’ .GBIOTB 

WRITE(2,*) 


MAXIMUM  NUMBER  OF  TIME  STEPS  =  ’ ,MAXT 
TOLERANCE  FOR  CONVERGENCE  =  ’ ,TOL 


2)  THEN 

TIME  STEP  ASSIGNED  TO  THE  "START  UP"’ 
PROBLEM  (  NO  EXTERNAL  DISTURBANCE  )  : 

BIOT  NUMBER  AT  THE  TOP  SURFACE  = 
BIOT  NUMBER  AT  THE  BOTTOM  SURFACE 


’  .GBIOT 
=  ’ .GBIOTB 


RATIO  (RHO*CP)GAGE  /  (RHO*CP)MAT.B 
RATIO  (K)KAT.B  /  (K)GAGE  *  ’ ,KBG 


' ,RCGB 


WRITE(2,*) ’ 
WRITER,*)’ 

ENDIF 

IF  ((TIME1  .EQ. 
WRITE(2,*) 
WRITE(2,*) ’ 
WRITE(2 ,*) ’ 
WRITE(2 ,*) ’ 
WRITE(2,*) ’ 
WRITE(2,*) ’ 
WRITE(2,*) ’ 
WRITE(2,*) ’ 


RATIO  (RHO*CP)GAGE  /  (RHO*CP)MAT.B  =  ’ , RCGB 
RATIO  (K)MAT.B  /  (K)GAGE  =  ’ ,KBG 

0)  .AND .  (MAXT  .GT.  D)  THEN 

TIME  STEPS  ASSIGNED  TO  THE  PROBLEM’ 

WITH  AN  EXTERNAL  DISTURBANCE  :  1  -  ’  ,  MAXT 

BIOT  NUMBER  AT  THE  TOP  SURFACE  *  ’ .DBIOT 
BIOT  NUMBER  AT  THE  BOTTOM  SURFACE  =  ’ ,DBIOTB 
RATIO  (QGEN/DBIQT) / (TFFINAL  -  TI)  *  ’ .GENRAT 
RATIO  (TFBACK  -  TI)/ (TFFRONT  -  TI)  «  ’ .BFRAT 
RATIO  (DBIOT)/(GBIOT)  *  ’ .DGBIRAT 
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WRITE(2,*) 
WRITE(2,*)  ’ 
WRITE (2,*)’ 
ENOIF 

IF  ((TIME1  .EQ. 
WRITEC2,*) 
WRITE(2,*)’ 
WRITE (2,*)* 
WRITE(2,*) * 
WRITE (2,*)’ 
WRITE(2,*) ’ 
WRITE(2,*) ’ 
WRITEC2,*)’ 
WRITE(2,*) 
WRITE(2 ,*) ’ 
WRITE(2,*) ’ 
ENDIF 

IF  ((TIME1  .GT. 
WRITE(2,«0 
WRITEC2.*)' 
WRITE(2 ,  *) * 
WRITEC2.*)’ 
WRITEC2 ,  *) ’ 
WRITEC2,*)' 
WRITE(2,*)’ 
WRITEC2,*)’ 
WRITE(2,*) 
WR1TE(2.*)’ 
WRITE(2,*)’ 
ENDIF 

IF  ((TIME1  .GT. 
WRITE(2,*) 
WRITEC2,*)’ 
WRITE(2,*) ’ 


RATIO  (RHO*CP)GAGE  /  (RHO*CP)MAT.B  »  ' , RCGB 
RATIO  (K)MAT.B  /  (K)GAGE  «  * ,KBG 

0)  .AND.  (NAXT  .EQ.  D)  THEN 

TIME  STEP  ASSIGNED  TO  THE  PROBLEM* 

WITH  AN  EXTERNAL  DISTURBANCE  :  1* 

BIOT  NUMBER  AT  THE  TOP  SURFACE  «  * .DBIOT 
BIOT  NUMBER  AT  THE  BOTTOM  SURFACE  *  *,DBIOTB 
RATIO  (QGEN/DBIOT)/(TFFINAL  -  TI) 

RATIO  (TFBACK  -  TI)/(TFFRONT  -  TI) 

RATIO  (DBIQT)/(GBIQT)  *  * .DGBIRAT 


’  .GENRAT 
*,BFRAT 


RATIO  (RHO*CP)GAGE  /  (RHO*CP)MAT.B  *  ’ ,RCGB 
RATIO  (K)MAT.B  /  (K)GAGE  =  \KBG 

0) .AND. (TIME1  .EQ.  MAXT})  THEN 

TIME  STEP  ASSIGNED  TO  THE  PROBLEM’ 

WITH  AN  EXTERNAL  DISTURBANCE  :  ’ .TIME1 

BIOT  NUMBER  AT  THE  TOP  SURFACE  =  ’ , DBIOT 
BIOT  NUMBER  AT  THE  BOTTOM  SURFACE  =  ’ .DBIOTB 
RATIO  (QGEN/DBIOT)/(TFFINAL  -  TI)  *  * .GENRAT 
RATIO  (TFBACK  -  TI)/ (TFFRONT  -  TI)  *  * , BFRAT 
RATIO  (DBIOT)/ (GBIOT)  =  * .DGBIRAT 


RATIO  (RHO*CP)GAGE  /  (RHO*CP)MAT.B 
RATIO  (K)MAT.B  /  (K)GAGE  *  ’ ,KBG 

0). AND. (TIME 1  .LT.  MAXT))  THEN 

TIME  STEPS  ASSIGNED  TO  THE  PROBLEM* 

WITH  AN  EXTERNAL  DISTURBANCE  :  ’ , 


’ .RCGB 


WRITE(2,*) 

WRITE(2,*) 

WRITE(2,*) 


TIME1 ,  ’  -  \MAXT 

BIOT  NUMBER  AT  THE  TOP  SURFACE 


.DBIOT 


BIOT  NUMBER  AT  THE  BOTTOM  SURFACE  *  DBIOTB 
RATIO  (QGEN/DBIOT)/(TFFINAL  -  TI)  =  ’.GENRAT 


WRITE(2,*) ’ 
WRITE(2,*) ’ 
WRITE(2 ,*) 
WRITE(2,*) ’ 
WRITE(2.*)’ 
ENDIF 


RATIO  (TFBACK  -  TI)/ (TFFRONT  -  TI)  =  \ BFRAT 
RATIO  (DBIOT)/ (GBIOT)  =  ’.DGBIRAT 

RATIO  (RHO*CP)GAGE  /  (RHO*CP)MAT.B  *  ’.RCGB 
RATIO  (K)MAT.B  /  (K)GAGE  «  ’ ,KBG 
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WRITE(2,*) 

HRITEC2.*) 

HRITEC2,*) ’RESULTS:' 

WRITE(2,*) 

C 

C . 

C  INITIALIZE  ARRAYS  AND  OTHER  PARAMETERS. 

C 

READ(3,*) 

DO  100  K  *  O.KMAX 
DO  110  I  *  O.IMAX 
THANEW(I.K)  *  0. 

READ (3, 5)  THAOLD(I.K) 

IF  (TIME1  .EQ.  0)  THAOLD(I.K)  *  THAOLD(I ,K)*GENRAT*DGBIRAT 
110  CONTINUE 

100  CONTINUE 
C 

DC  120  I  =  0, IGAGE 

IF  (I  .LE.  IGENMX)  THEN 
I GEN (I)  =  1. 

ELSE 

IGEN(I)  =  0. 

END  IF 
120  CONTINUE 
C 

DELR  =  1 . /FLOAT(IGAGE) 

DELZ  «  1 ./FLOAT(KMAX) 

C 

LAMR1  «  DELT/DELR**2 
LAMZ1  »  DELT/DELZ**2 
LAMZ2  *  DELT/DELZ 
C 

NUMER  *  2 . *RCGB*FLOAT ( IGAGE) 

DENOM  =  RCGB* (FLOAT (IGAGE )  -  .25)  ♦  (FLOAT( IGAGE)  ♦  .25) 

RCGGB  *  NUMER/DENOM 
C 

NUMER  »  (FLOAT (IGAGE)  -  .25)  ♦  KBG* (FLOAT (IGAGE)  ♦  .25) 

DENOM  =  2 . *FLOAT(IGAGE) 

KGBG  *  NUMER/DENOM 
C 

PI  »  (FLOAT(IMAX)  -  .5)/(FL0AT(IMAX)  -  .25) 

P2  ■  1.  -  l./(2.*FL0AT(IGAGE)) 

P3  ■  1.  ♦  l./(2.*FL0AT(IGAGE)) 

C 

C . - . - . - . . . 
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C  I************************************************************* 

c . 

C  MAIN  LOOP. 

C 

IF  (TIME1  .EQ.  0)  THEN 
BIOT  ■  DBIOT 
BIOTB  «  DBIOTB 
ELSE 

BIOT  ■  GBIOT 
BIOTB  «  GBIOTB 
ENDIF 
COUNT  =  0 
MAXCH  «  99.99 
MAXDER  *  0. 

C 

C . - . . 

DO  SOO  N  =  l.MAXT 
IF  ((MAXCH/DELT)  .GT.  TOL)  THEN 

C . . . . . 

c 

c 

c  . - . . 

C  REDEFINE  BIOT  NUMBERS  AND  UPDATE  THAOLD  WHEN  STARTING 

C  THE  DISTURBANCE  PROBLEM. 

C 

IF  (N  .EQ.  TIMED  THEN 
BIOT  «  DBIOT 
BIOTB  =  DBIOTB 
DO  1000  I  *  O.IMAX 

DO  1010  K  *  O.KMAX 

THAOLD(I.K)  =  THAOLDCI ,K)*GENRAT*DGBIRAT 
1010  CONTINUE 

1000  CONTINUE 

ENDIF 
C 

C  . - . — . 

C  SOLVE  FOR  NEXT  TIME  STEP. 

C 

C 

C  FRONT  OF  GAGE 

C 
C 
C 

C  (A) 

C 
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A  •  4 . * (LRRAT**2) *LAMRi*  THAOLD(l.O) 

B  ■  2.*LAMZ1*  THA0LD(0,1) 

Cl  «  1.  -  4 . * (LRRAT**2) *LAMR1  -  2.*LAMZ1  -  2.*LAMZ2*BI0T 
C  «  Cl*  THAOLD(O.O) 

IF  (N  .LT.  TIMED  THEN 

D  «  2 . *LAMZ2*BI0T*  IGEN(O) 

ELSE 

D  »  2.*LAMZ2*BI0T*(1 .  +  IGEN(0)*GENRAT) 

ENDIF 

THANEW(O.O)  “A+B+C+D 
C 

C  CB) 

C 

DO  640  I  *  l.IGAGE-l 

A1  *  1.  -  l./(2.*FL0AT(I)) 

A2  *  1.  ♦  1 . / (2 . *FLOAT(I) ) 

A3  =  (LRRAT**2) *LAMR1 

A  *  A3*(  Al*  THA0LD(I-1 ,0)  ♦  A2*  THAOLD(I-H.O)  ) 

B  *  2.*LAMZ1*  THAOLD(I.l) 

Cl  *  1.  -  2 . * (LRRAT**2) *LAMR1  -  2.*LAMZ1  -  2.*LAMZ2*BI0T 
C  *  Cl*  THAOLD(I.O) 

IF  (N  .LT.  TIMED  THEN 

D  *  2.*LAMZ2*BI0T*  IGEN(I) 

ELSE 

D  *  2 . *LAMZ2*BI0T*( 1 .  +  I GEN ( I ) *GENRAT ) 

ENDIF 

THANEW(I.O)  =A*B+C+D 
640  CONTINUE 

C 

C  (C) 

C 

Al  *  RCGGB*(LRRAT**2)*LAMR1 

A  *  Al*(  P2*  THA0LD(IGAGE-1 ,0)  ♦  P3*KBG*  THA0LD(IGAGE+1 ,0)  ) 
B  *  2 . *RCGGB*KGBG*LAMZ 1 *  THAOLDQGAGE,  1) 

Cl  *  -  RCGGB* (LRRAT**2) *P2*LAMR1 

C2  «  -  RCGGB*KBG*(LRRAT**2)*P3*LAMR1 

C3  *  -  2 . *RCGGB*KGBG*LAMZ1  -  2.*RCGGB*BI0T*LAMZ2 

C  «  (1.  *  Cl  ♦  C2  ♦  C3)*  THAOLD(IGAGE.O) 

IF  (N  .LT.  TIMED  THEN 

D  =  2 . *RCGGB*BI0T*LAMZ2*  IGEN(IGAGE) 

ELSE 

D  =  2.*RCGGB*BI0T*LAMZ2*(1 .  +  I GEN ( IGAGE) *GENRAT) 

ENDIF 

THANEW(IGAGE.O)  »A+B+C*D 
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c 

c 

c 


CD) 


DO  660  I  -  IGAGE+1 , IMAX-1 
A1  »  1.  -  l./(2.*FL0AT(I)) 

A2  »  1.  ♦  l./(2.*FL0AT(I)) 

A3  ■  RCGB*KBG*(LRRAT**2)*LAMR1 
A  «  A3*(  Al*  THAOLD(I-l.O)  ♦  A2*  THAOLDU+l.O)  ) 

B  ■  2 . *RCGB*KBG*LAMZ1*  THAOLD(I.l) 

Cl  *  -  2 . *RCGB*KBG* (LRRAT**2) *LAMR1 

C2  *  -  2 . *RCGB*KBG*LAMZ 1 

C3  »  -  2 . *RCGB*BI0T*LAMZ2 

C  *  (1.  ♦  Cl  +  C2  +  C3)*  THAOLD(I.O) 

IF  (N  .LT.  TIMED  THEN 
D  =  0. 

ELSE 

D  *  2 . *RCGB*BI0T*LAMZ2 
ENDIF 

THANEWCI ,0)  =A+B*C+D 
660  CONTINUE 

C 

C  (E) 

C 

A  =  2 . *RCGB*KBG* (LRRAT**2) *P1*LAMR1*  THA0LD(IMAX-1 ,0) 
B  *  2 . *RCGB*KBG*LAMZ1*  THAOLD(IMAX.l) 

Cl  *  -  2.*RCGB*KBG*(LRRAT**2)*P1*LAMR1 
C2  »  -  2.*RCGB*KBG*LAHZ1 
C3  *  -  2 . *RCGB*BI0T*LAMZ2 
C  ■  (1.  ♦  Cl  ♦  C2  ♦  C3)*  THAOLD(IMAX.O) 

IF  CN  .LT.  TIMED  THEN 
D  «  0. 

ELSE 

D  *  2 . *RCGB*BI0T*LAMZ2 
ENDIF 

THANEW(IMAX.O)  -A  +  B  +  C  +  D 
C 
C 

C  INTERIOR  OF  GAGE 

C 

C 

DO  700  K  «  l.KMAX-1 
(F) 


C 

C 

C 
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A  ■  4.*(LRRAT**2)*LAMR1*  THAOLD(l.K) 

B  -  LAMZ1*(  THAOLD(0,K-1)  ♦  THAOLD(0,K-H)  ) 

Cl  -  1.  -  4 . * (LRRAT**2) *LAMR1  -  2.*LAMZ1 
C  »  Cl*  THAOLD(O.K) 

THANEV(O.K)  ■  A  ♦  B  ♦  C 
C 

C  (G) 

C 

DO  740  I  *  l.IGAGE-1 

A1  *  1.  -  l./(2.*FL0AT(I)) 

A2  *  1.  ♦  l./(2.*FL0AT(I)) 

A3  *  (LRRAT**2) *LANR1 

A  »  A3*(  Al*  THAOLD(I-l.K)  ♦  A2*  THA0LD(I+1 ,K)  ) 

B  =  LAMZ1*(  THAOLD(I ,K-1)  ♦  THAOLDd.K+l)  ) 

Cl  =  1.  -  2 . * (LRRA1 **2) *LAMR1  -  2.*LAMZ1 
C  =  Cl*  THAOLD(I.K) 

THANEW(I.K)  -  A  ♦  B  ♦  C 
740  CONTINUE 

C 

C  (H) 

C 

Al  *  RCGGB*(LRRAT**2)*LAMR1 

A  -  Al*(  P2*  THAOLD  dG  AGE- 1,K)  +  P3*KBG*  THAOLD  dGAGE+ 1  ,K)  ) 
B1  *  RCGGB*KGBG*LAMZ1 

B  a  Bl*(  THAOLD ( IGAGE , K- 1 )  +  THAOLD ( I GAGE ,  K*  1 )  ) 

Cl  »  -  RCGGB* (LRRAT**2) *P2*LAMR1 
C2  *  -  RCGGB*KBG* (LRRAT**2) *P3*LAHR1 
C3  *  -  2 . *RCGGB*KGBG*LAMZ1 
C  *  (i.  ♦  Cl  ♦  C2  ♦  C3)*  THAOLD ( I GAGE, K) 

THANEW( IGAGE, K)  -  A  +  B  ♦  C 
C 

C  (I) 

C 

DO  760  I  =  IGAGE+1 , IMAX-1 
Al  *  1.  -  l./(2.*FL0AT(I)) 

A2  «  1.  +  l./(2.*FL0AT(I)) 

A3  a  RCGB*KBG*(LRRAT**2) *LAMR1 
A  «  A3*(  Al*  THAOLD(I-l.K)  +  A2*  THAOLDd+1 ,K)  ) 

B1  =  RCGB*KBG*LAMZ1 

B  «  Bl*(  THAOLD(I.K-l)  ♦  THAOLD(I ,K+l)  ) 

Cl  =  -  2 . *RCGB*KBG* (LRRAT* *2) *LAMR1 
C2  «  -  2 . *RCGB*KBG*LAMZ1 
C  *  (1.  ♦  Cl  ♦  C2)*  THAOLD(I ,K) 

THANEWd  ,K)  ■  A  +  B  ♦  C 
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760  CONTINUE 

C 

C  (J) 

C 

A  *  2 . *RCGB*KBG* (LRRAT**2) *P1*LAMR1*  THA0LD(IMAX-1,K) 

B1  •  RCGB*KBG*LAMZ1 

B  ■  8l*(  THA0LD(IMAX,K-1)  4  THAOLD ( IMAX ,  K+ 1)  ) 

Cl  «  -  :>.*RCGB*KBG*(LRRAT**2)*P1*LAMR1 
C2  *  -  2 . *RCGB*KBG*LAMZ 1 
C  ■  (1.  ♦  Cl  ♦  C2)*  THAOLD ( IMAX, K) 

THANEW(IMAX.K)  *  A  ♦  B  ♦  C 
C 

700  CONTINUE 

C 

C 

C  BACK  OF  GAGE 

C 
C 
C 

C  (K) 

C 

A  =  4 . * (LRRAT**2) *LAMR1*  THA0LD(1 ,KMAX) 

B  *  2 .*LAMZ1*  THAOLD (O.KMAX-1) 

Cl  *  1.  -  4 . * (LRRAT**2) *LAMR1  -  2.*LAMZl  -  2 . *LAMZ2*BI0TB 
C  =  Cl*  THAOLD(O.KMAX) 

IF  (N  .LT.  TIMED  THEN 
D  «  0. 

ELSE 

D  *  2 . *LAMZ2*BI0TB*  BFRAT 
END  IF 

THANEW(O.KMAX)  =A*B+C*D 
C 

C  (L) 

C 

DO  840  I  •  1 .IGAGE-l 

A1  =  1.  -  l./(2.*FL0AT(I)) 

A2  *  1.  *  l./(2.*FL0AT(I)) 

A3  =  (LRRAT**2) *LAMR1 

A  =  A3*(  Al*  THAOLD(I-l.KMAJ)  4  A2*  THA0LD(I+1 ,KMAX)  ) 

B  =  2.*LAMZ1*  THAOLD(I .KMAX-l) 

Cl  *  1.  -  2 . * (LRRAT**2) *LAMR1  -  2.*LAMZ1  -  2 . *LAMZ2*BI0TB 
C  «  Cl*  THAOLD ( I, KMAX) 

IF  (N  .LT.  TIMED  THEN 
D  -  0. 
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ELSE 

D  •  2 .  *LAMZ2*BI0TB*  BFRAT 
ENDIF 

THANEV(I.KMAX)  -A  +  B  +  C  +  D 
840  CONTINUE 

C 

C  (M) 

C 

A1  ■  RCGGB*  (LRRAT**2)  *P2*LAMR1 
A2  «  RCGGB*KBG*(LRRAT**2)*P3*LAMR1 

A  *  Al*  THA0LD(IGAGE-1 ,KMAX)  ♦  A2*  THA0LD(IGAGE+1 ,KMAX) 

B  *  2.*RCGGB*KGBG*LAMZ1*  THAOLD(IGAGE.KMAX-l) 

Cl  *  -  RCGGB* (LRRAT**2)*P2*LAMR1 

C2  *  -  RCGGB*KBG*(LRRAT**2)*P3*LAMR1 

C3  =  -  2 . *RCGGB*KGBG*LAMZ 1  -  2.*RCGGB*BI0TB*LAMZ2 

C  *  (1.  ♦  Cl  ♦  C2  +  C3)*  THAOLD(IGAGE.KMAX) 

IF  (N  .LT.  TIMED  THEN 
D  =  0. 

ELSE 

D  =  2 .  *RCGGB*BI0TB*LAMZ2*  BFRAT 
ENDIF 

THANEW(IGAGE.KMAX)  -A+B+C+D 
C 

C  (N) 

C 

DO  860  I  =  IGAGE+1 ,IMAX-1 
Al  =  1.  -  l./(2.*FL0AT(I)) 

A2  =  1.  ♦  l./(2.*FL0AT(I)) 

A3  =  RCGB*KBG*(LRRAT**2)*LAMR1 

A  =  A3*(  Al*  THA0LD(I-1,KHAX)  ♦  A2*  TH AOLD ( I ♦ 1 , KMAX )  ) 
B  *  2.*RCGB*KBG*LAMZ1*  THAOLD(I .KMAX-l) 

Cl  *  -  2 . *RCGB*KBG* (LRRAT**2) *LAMR1 
C2  *  -  2 . *RCGB*KBG*LAMZ 1 
C3  =  -  2 . *RCGB*BI0TB*LAMZ2 
C  *  (1.  ♦  Cl  ♦  C2  ♦  C3)*  THAOLD(I.KMAX) 

IF  (N  .LT.  TIMED  THEN 
D  =  0. 

ELSE 

D  =  2 . *RCGB*BI0TB*LAMZ2*  BFRAT 
ENDIF 

THANEW(I.KMAX)  -A+B+C+D 
860  CONTINUE 

C 

C  (0) 
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c 

A  ■  2 . *RCGB*KBG* (LRRAT**2) *P1*LAMR1*  THA0LD(IMAX-1 , KMAX) 

B  «  2 . *RCGB*KBG*LAMZ1*  THAOLD (IMAX. KMAX- 1) 

Cl  -  -  2 . *RCGB*KBG* (LRRAT**2) *P1*LAMR1 
C2  «  -  2 . *RCGB*KBG*LAMZ 1 
C3  *  -  2 . *RCGB*BI0TB*LAMZ2 
C  ■  (1.  ♦  Cl  ♦  C2  ♦  C3)*  THAOLD(IMAX.KMAX) 

IF  (N  .LT.  TIMED  THEM 
D  «  0. 

ELSE 

D  «  2 . *RCGB*BI0TB*LAMZ2*  BFRAT 
END  IF 

THANEW(IMAX.KMAX)  =  A  +  B  +  C*D 
C 
C 
C 

C  . . . . 

C  WRITE  NON-DIMENSIONAL  TEMPERATURE  VALUES  AT  THE  HEATED  DISK 

C  TO  FILE  "DISKT"  (DURING  DISTURBANCE  PROBLEM  ONLY). 

C  **  MODIFIED  TO  WRITE  TO  "DISKT"  ONLY  AT  INITIAL  TIME  STEP 

C  OF  DISTURBANCE  PROBLEM  AND  AT  TIME  STEPS  T+  =  .001-. 070 

C  WITH  DELT=.001  AFTER  THE  DISTURBANCE  PROBLEM  STARTS. 

C  BE  CAREFUL  TO  HAVE  ENOUGH  TIME  STEPS  TO  HAVE  ATLEAST 

C  .070  IN  NON-DIMENSIONAL  TIME  FOR  THE  DISTURBANCE  (I.E. 

C  NEED  TO  HAVE  COUNT  ATLEAST  EQUAL  TO  70  WHEN  THE  PROGRAM 

C  FINISHES) . 

C 

IF  (N  .GE.  TIMED  THEN 

IF  (  (N.Eq. TIMED  .OR.  ( (TIME1 . EQ . 0) . AND . (N . EQ . 1 ) )  )  THEN 
COUNT  *  0 

WRITE (10. 25)  COUNT 
DO  1100  I  =  0,1 GEN MX 

WRITE(10,30)  I.THAOLD(I.O) 

1100  CONTINUE 

END  IF 
C 

X  *  N  -  MAX(TIMEl.l)  +  1 
Y  =  N INT ( . 00 1 /DELT ) 

IF  (  MOD(X.Y)  .EQ.  0  )  THEN 
COUNT  *  COUNT  ♦  1 
WRITE(10,25)  COUNT 
DO  1125  I  -  O.IGENMX 

WRITE( 10,30)  I.THANEW(I.O) 
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1125  CONTINUE 

ENDIF 
END  IF 
C 

C  . 

C  ROLL  DOWN  THAOLD(l :20, 1 :20)  AND  FIND  MAXIMUM  VISUAL 

C  PERCENT  CHANGE.  (BE  CAREFUL  NOT  TO  DIVIDE  BY  ZERO.) 

C 

MAXTHA  -  0. 

DO  900  I  *  O.IMAX 
DO  910  K  «  O.KMAX 

IF  (  ABS (THANEW ( I ,  K  )  )  .GT.  MAXTHA)  THEN 
MAXTHA  «  ABS (THANEW ( I, K)) 

ENDIF 

910  CONTINUE 

900  CONTINUE 

C 

MAXCH  =  0. 

DO  950  I  *  O.IMAX 
DO  960  K  =  O.KMAX 

CHNG  =  ABS(  THANEW(I.K)  -  THAOLD(I.K)  ) 

IF  (  MAXTHA  .GT.  l.E-15)  THEN 
VISCHNG  *  CHNG/MAXTHA 
IF  (VISCHNG  .GT.  MAXCH)  MAXCH  *  VISCHNG 
ELSE 

IF  (CHNG  .GT.  MAXCH)  MAXCH  =  CHNG 
ENDIF 

THAOLD(I.K)  *  THANEW(I.K) 

960  CONTINUE 

950  CONTINUE 

C 

C  . . . 

C  FIND  MAXIMUM  RADIAL  DIRECTION  DERIVATIVE  AT  OUTER 

C  BOUNDARY  OF  THE  PROBLEM.  SAVE  THE  MAXIMUM  VALUE 

C  FROM  THE  WHOLE  RUN  (ALL  ITERATIONS  INCLUDED)  TO 

C  PRINT  OUT  WHEN  FINISHED  WITH  THE  MAIN  LOOP. 

C 

DO  1200  K  =  O.KMAX 

DER  »  (THANEW(IMAX.K)  -  THANEW (IMAX- 1 ,K))/DELR 
IF  (ABS (DER)  .GT.  ABS(MAXDER) )  MAXDER  =  DER 

1200  CONTINUE 

C 

C . . 

ELSE 
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C . - . 

c 

C  CONVERGES.  JUMP  OUT  OF  LOOP. 

C 

WRITE(2,*) 

WRITE (2,*) ‘CONVERGES  AFTER  »,N-1,‘  ITERATIONS.' 

WRITE(2f*) ’MAXIMUM  CHANGE  ON  THE  LAST  ITERATION  WAS  \MAXCH,’.' 
GOTO  2000 
C 

C . 

ENDIF 

500  CONTINUE 

C . 

c 

C  IF  YOU  EXECUTE  THESE  STATEMENTS,  THEN  YOU  COMPLETED  THE 

C  MAXIMUM  NUMBER  OF  ITERATIONS  FOR  THE  LOOP  WITHOUT  CONVERGING. 

C  PRINT  MESSAGE  IF  EXPECTING  SOME  CONVERGENCE  (I.E.  TOL  >  0). 

C 

IF  (TOL  .GT.  0.)  THEN 
WRITE(2,*) 

WRITE (2,*) 'QUITS  WITHOUT  CONVERGING  AFTER  \N-1,’  ITERATIONS.' 
WRITE(2,*) 'MAXIMUM  CHANGE  ON  THE  LAST  ITERATION  WAS  ’ ,MAXCH ,  ’ .  ’ 
ENDIF 
C 

C . 

c  ******************************************************************* 

C . 

2000  CONTINUE 
C 

C . - . - . - . 

C  PRINT  VALUE  OF  RADIAL  DIRECTION  DERIVATIVE  AT  OUTER  BOUNDARY  THAT 
C  WAS  THE  LARGEST  (IN  ABSOLUTE  VALUE)  DURING  THE  RUN. 

C 

WRITE(2,«0 

WRITE(2,*) 'MAXIMUM  RADIAL  DIRECTION  DERIVATIVE  AT  OUTER  BOUNDARY’ 
WRITE(2,*) 'DURING  THE  RUN  WAS  ’ .MAXDER, ’ . ’ 

C 

C . - . - . 

C  PRINT  CHECK  ON  VALUE  OF  COUNT. 

C 

IF  (COUNT  .LT.  70)  THEN 
WRITE (2,*) 

WRITE(2,*) 'WARNING:  COUNT  LESS  THAN  70.  COUNT  =  ’.COUNT 
WRITE(2,*) 

ELSE 
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WRITE(2,*) 


WRITE(2,*) ’VALUE  OF  COUNT  IS  * .COUNT 
WRITE (2,*) 

WRITE(2,*) 

ENDIF 

C 

. . 

C  FORMAT  STATEMENTS. 

C 

5  FORMAT (20X ,E17 . 10) 

25  FORMAT (IS) 

30  F0RMAT(I5,SX,E17.10) 

C 

END 
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D.S  Code  to  Compare  the  Finite- Difference  and  the  Series  Solution  Estimates  for 
Surface  Heat  Flux 
PROGRAM  HTTR2H 

C  **************************************************************** 


C  *  2LT  JOSEPH  A.  BONAFEDE,  GA-88M  * 

C  *  FALL  1987  * 

C  *  * 

C  *  ADVISOR:  DR.  JAMES  E.  HITCHCOCK  * 


C  **************************************************************** 

c 

C  DRIVER  FOR  SUBROUTINE  GETQ. 

C  GETS  INPUTS  FOR  CALL  TO  SUBROUTINE  GETQ  FROM  THE  INPUT  FILE  TO 
C  GAGE2,  I.E.  "G2INP"  .  OPENS  THE  OUTPUT  FILE  FOR  GETQ  WHICH  IS 
C  ALSO  THE  OUTPUT  FILE  FOR  GAGE2,  I.E.  "G20UT"  (UNIT=2) . 

C  OPENS  THE  INPUT  FILE  FOR  GETQ,  I.E.  "DISKT"  (UNIT=10) . 

C  OPENS  THE  OUTPUT  DATA  FILES  "PTVDF"  (UNIT=11),  "PTPDF"  (UNIT=12) , 
C  "PTSER"  (UNIT=13) ,  AND  "PTFD"  (UNIT=14) . 

C 

C  **  MODIFIED  TO  WORK  WITH  PROGRAMS  GAGE2M  OR  GAGE2N  AND 
C  **  SENDS  VALUE  OF  70  FOR  COUNT  AND  .001  FOR  DELT  ALWAYS. 

C 

C . 

C  DECLARE  VARIABLES. 

C 

IMPLICIT  CHARACTER(A-Z) 

INTEGER  IGENMX, COUNT 
REAL  DB I OT , DELT , GENRAT 
INTEGER  MAXT.TDIST 
C 

C . - . - . . . 

C  OPEN  INPUT  AND  OUTPUT  FILES. 

C 

OPEN (UNIT® 1 , FILE® ’ G2INP ’ , STATUS® » OLD » ) 

OPEN (UNIT®2 ,  FILE® ' G2QUT * , STATUS® » OLD  * ) 

OPEN (UNIT® 10 , FILE® ’ DISKT ’ , STATUS* ’ OLD ’ ) 

OPEN (UNIT* 1 1 .FILE® » PTVDF ’ , STATUS* * NEW • ) 

OPEN (UNIT® 1 2 , F I LE= ’ PTPDF ’ , STATUS*  *  NEW » ) 

OPEN (UNIT® 13 , FILE® ’ PTSER ’ , STATUS* ’ NEW * ) 

OPEN (UN IT® 14 , F I LE= * PTFD STATUS® ’ NEW ’ ) 

REWIND (UNIT® 1) 

C 

C . . . . . 

C  GET  INPUTS  TO  SUBROUTINE. 

C 
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READCl,*) 
READCl,*) 
READCl,*)  MAXT 
READCl,*)  TDIST 
READCl,*) 
READCl,*) 
READCl,*) 
READCl,*)  IGENMX 


READCl,*) 

READCl,*) 

READCl,*) 

READCl,*)  DBIOT 
READCl,*) 

READCl.*)  GENRAT 
C 

IF  C TDIST  .GT.  MAXT)  THEN 
WRITEC2,*) 

WRITEC2,*) 'PROBLEM:  TDIST  GREATER  THAN  MAXT.  THIS  WAS  NOT’ 
WRITEC2,*) ’A  DISTURBANCE  PROBLEM.’ 

STOP 

ENDIF 

C 

DELT  «  .001 
COUNT  =  70 
C 

C . . . . . . 

C  CALL  SUBROUTINE. 

C 

CALL  GETQ  C  COUNT, IGENMX, DBIOT, DELT, GENRAT  ) 

C 

END 

C  ****************************************************************** 

c  ****************************************************************** 
SUBROUTINE  GETQ  C  NUMT.IMX, BIOT, DELT, GENRAT  ) 

C 

c  ***************************************************************** 


C  *  2LT  JOSEPH  A.  BONAFEDE,  GA-88M  * 

C  *  FALL  1987  * 

C  *  * 

C  *  ADVISOR:  DR.  JAMES  E.  HITCHCOCK  * 


C  ***************************************************************** 

c 

C  THIS  SUBROUTINE  DETERMINES  THE  NON-DIMENSIONAL  HEAT 
C  TRANSFER  RATE  AT  THE  HEATED  DISK  FOR  BOTH  THE  SERIES 
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C  SOLUTION  AND  THE  FINITE  DIFFERENCE  SOLUTION. 

C  THE  SUBROUTINE  COMPARES  THE  DIFFERENT  VALUES  AND  PRINTS 
C  THE  RESULTS  IN  THE  FILE  "G20UTn  (UNIT-2). 

C  DATA  NEEDED  TO  DETERMINE  THE  NON-DIMENSIONAL  HEAT  TRANSFER 
C  RATES  IS  STORED  IN  THE  FILE  "DISKT"  (UNIT-10). 

C 

C  GLOSSARY  OF  MAIN  VARIABLES: 

C  NUMT  -  THE  NUMBER  OF  SURFACE  HEAT  TRANSFER  VALUES  TO  CALCULATE. 

C  INPUT 

C  IMX  -  THE  MAXIMUM  VALUE  FOR  RADIAL  NODE  INCLUDED  IN  THE  DISK, 

C  I.E.  IGENMX.  (NODES  I-O, IMX  ARE  INCLUDED  IN  THE  DISK.) 

C  INPUT 

C  BIOT  -  THE  SURFACE  BIOT  NUMBER  FOR  THE  DISTURBANCE  PROBLEM. 

C  INPUT 

C  DELT  -  THE  NON-DIMENSIONAL  TIME  STEP  USED  IN  THE  PROGRAM. 

C  INPUT 

C  GENRAT  -  THE  NON-DIMENSIONAL  RATIO  (QG/DBIOT)/(TFFINAL  -  TI) 

C  USED  IN  THE  PROGRAM 

C  INPUT 

C  THANOH  -  MATRIX  STORING  THE  VALUES  FOR  THATA(I)  AT  TIME  =  J 
C  FOR  I-O. IMX. 

C  USED  INTERNALLY 

C  THAPRE  -  MATRIX  STORING  THE  VALUES  FOR  THATA(I)  AT  TIME  T  *  J-l 
C  FOR  1=0, IMX. 

C  USED  INTERNALLY 

C  SUM  -  MATRIX  WHICH  STORES  THE  VALUES  (THANOW(I)-THAPRE(I))/B 
C  FOR  I-O, IMX.  B  IS  A  FUNTION  OF  TIME  T  (TIME  OF  THE  SURFACE 

C  HEAT  TRANSFER  VALUE)  AND  TIME  J  (TIME  ASSOCIATED  WITH  THANOW) . 

C  USED  INTERNALLY 

C  B  -  THE  DENOMINATOR  TERM  INSIDE  THE  SERIES  FOR  THE  SERIES  SOLUTION. 

C  B  IS  A  FUNCTION  OF  TIME  T  (TIME  OF  THE  SURFACE  HEAT  TRANSFER 

C  VALUE)  AND  TIME  J  (TIME  ASSOCIATED  WITH  THANOW) . 

C  USED  INTERNALLY 

C  CONS  -  THE  CONSTANT  TERM  IN  THE  SERIES  SOLUTION. 

C  USED  INTERNALLY 

C  CONFD  -  THE  CONSTANT  TERM  IN  THE  FINITE  DIFFERENCE  SOLUTION. 

C  USED  INTERNALLY 

C  XS  -  VARIABLE  USED  FOR  AN  INTERMEDIATE  RESULT  IN  THE  SERIES  SOLUTION. 
C  USED  INTERNALLY 

C  XFD  -  VARIABLE  USED  FOR  AN  INTERMEDIATE  RESULT  IN  THE  FINITE 

C  DIFFERENCE  SOLUTION. 

C  USED  INTERNALLY 

C  QS  -  THE  FINAL  VALUE  FOR  AVERAGE  SURFACE  HEAT  TRANSFER  AT  THE  DISK 
C  FOR  THE  SERIES  SOLUTION. 
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C  USED  INTERNALLY 

C  QFD  -  THE  FINAL  VALUE  FOR  AVERAGE  SURFACE  HEAT  TRANSFER  AT  THE  DISK 
C  FOR  THE  FINITE  DIFFERENCE  SOLUTION. 

C  USED  INTERNALLY 

C  VISDIFF  -  THE  PERCENT  OF  THE  FULL  SCALE  VALUE  (MAX  THEORETICAL  VALUE 
C  FOR  SURFACE  HEAT  TRANSFER)  BY  WHICH  THE  RESULTS  FOR 

C  AVERAGE  SURFACE  HEAT  TRANSFER  AT  THE  DISK  DIFFER. 

C  (FULL  SCALE  -  1.  ♦  GENRAT) . 

C  USED  INTERNALLY 

C  PERDIFF  -  THE  PERCENT  DIFFERENCE  BETWEEN  THE  RESULTS  FOR  AVERAGE 
C  SURFACE  HEAT  TRANSFER  AT  THE  DISK  (USING  THE  FINITE 

C  DIFFERENCE  VALUE  IN  THE  DENOMINATOR). 

C  USED  INTERNALLY 

C  T  -  INDEX  DENOTING  THE  TIME  STEP  FOR  WHICH  WE  ARE  CURRENTLY 
C  CALCULATING  THE  SURFACE  HEAT  TRANSFER. 

C  USED  INTERNALLY 

C  J  -  INDEX  DENOTING  THE  TIME  STEP  INSIDE  THE  SERIES  (ASSOCIATED 
C  WITH  THANOW) . 

C  USED  INTERNALLY 

Cl-  INDEX  DENOTING  THE  RADIAL  NODE  (1=0, IMX). 

C  USED  INTERNALLY 

C  PI  -  THE  CONSTANT  PI  (3.141  ...) 

C  USED  INTERNALLY 

C  DECLARE  VARIABLES: 

IMPLICIT  CHARACTER ( A- Z) 

INTEGER  NUMT.IMX 
REAL  BIOT, DELT, GENRAT 

REAL  SUM ( 1 : 100) , THAPRE (1:100), THANOW (1:100) 

REAL  PI.CONS.CONFD 
REAL  B,XS,XFD,QS,QFD 
REAL  VISDIFF, PERDIFF 
INTEGER  T,I, J 
C 

C . - . - . — 

C  GET  PARAMETERS . 

C 

PI  =  4. *ATAN(1 . ) 

CONS  =  2 ,/(BIOT*SQRT(DELT*PI)*( (FLOAT(IMX)  +  .5)**2)) 

CONFD  «  1 ./((FLOAT(IMX)  ♦  .5)**2) 

C 

C  — . - . - - - - — 

C  PRINT  HEADER  TO  "G20UT"  (UNIT  2) . 

C 

WRITE(2,50) 
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C . 

C  GET  RESULTS  FOR  EACH  TIME  T. 
C 


DO  500  T  -  l.NUMT 
C 

REWIND (UNIT* 10) 

C 

C  /*READ  VALUES  AT  TIME  ZERO.*/ 

C 

READ (10,*) 

DO  600  I  =  O.IMX 

READ(10,25)  THAPRE(I) 

600  CONTINUE 

C 

C  /*GET  SUM(I)  FROM  J=1,T-1  FOR  ALL  I*/ 

C 

DO  700  I  *  O.IMX 
SUM(I)  =  0. 

700  CONTINUE 

DO  800  J  =  l.T-l 

B  =  SQRT(FLOAT(T-J) )  +  SQRT(FLOAT(T- J+l) ) 
READ(10,*) 

DO  850  I  *  O.IMX 

READ (10, 25)  THANOW(I) 

SUM(I)  =  SUM(I)  ♦  (THANOW(I)  -  THAPRE(I))/B 


THAPRE(I)  =  THANOW(I) 

850  CONTINUE 

800  CONTINUE 

C 

C  /*GET  FINAL  VALUE  OF  SUM(I)  FOR  EACH  I  (INCLUDES  TERM  FOR*/ 

C  /*  J=T) .  ALSO,  GET  FINAL  VALUES  FOR  QS  AND  QFD .  */ 

C 

QS  *  0. 

QFD  *  0. 

READ(10,*) 

DO  900  I  =  O.IMX 

READ (10, 25)  THANOW(I) 

SUM(I)  =  SUM(I)  +  (THANOW(I)  -  THAPRE(I)) 

IF  (I  .EQ.  0)  THEN 
XS  =  ,25*SUM(I) 

XFD  *  . 25*(1 .  -  THANOW(I)) 

ELSE 

XS  =  2 . *FLOAT(I)*SUM(I) 

XFD  *  2.*FL0AT(I)*(1 .  -  THANOW(I)) 
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ENDIF 

QS  «  QS  ♦  XS 
QFD  *  QFD  ♦  XFD 
900  CONTINUE 

QS  «  QS  *  CONS 
QFD  -  QFD  *  CONFD 
C 

C  /^COMPARE  QS  AND  QFD  AND  PRINT  RESULTS  TO  ,,G20UT,,  (UNIT  2)*/ 

C  /*AND  "PTDIFF"  (UNIT  11)  */ 

C 

IF  (  ABS(QFD)  .GT.  IE-15  )  THEN 

VISDIFF  *  ((QS  -  QFD)*100. )/(l .  ♦  GENRAT) 

PERDIFF  =  ((QS  -  QFD) *100 . )/ABS(QFD) 

WRITE(2,60)  DELT*FLOAT(T) , QS.QFD, VISDIFF, PERDIFF 
WRITE(11 ,80)  DELT*FLOAT(T) .VISDIFF 
WRITE(12,80)  DELT*FLOAT(T) .PERDIFF 
WRITE(13,80)  DELT*FLOAT(T) ,QS 
WRITE(14,80)  DELT*FLOAT(T) ,QFD 
ELSE 

VISDIFF  =  ((QS  -  qFD)*100.)/(l.  +  GENRAT) 

WRITE(2 ,70)  DELT*FLOAT(T) .QS.QFD, VISDIFF 
WRITE(11 ,80)  DELT*FLOAT(T) .VISDIFF 
WRITE(12 ,85)  DELT*FLOAT(T) 

VRITEd3,80)  DELT*FLOAT(T)  ,QS 
WRITE (14, 80)  DELT*FLOAT(T) ,QFD 
ENDIF 
C 

500  CONTINUE 

C 

C . 

C  FORMAT  STATEMENTS. 

25  FORMAT (10X.E17 . 10) 

50  FORMAT ('NON-DIM.  TIME’ ,T21 , ’HEAT  TRANSFER’ ,T41 , ’HEAT  TRANSFER’/ 
ft  ’ (AFTER  DISTURBANCE) ’ ,T21 , ’ (SERIES  SOLUTION) ’ ,T41 , 

t  ’(FINITE  DIFF  SOLU) ’ ,T61 , ’VISUAL  DIFFERENCE’ ,T81 , 

t  ’PERCENT  DIFFERENCE’) 

60  F0RMAT(F10.7,10X,4(E17. 10, 3X) 

70  FORMAT(F10.7, 10X,3(E17 . 10, 3X) , ’  NOT  COMPARED’) 

80  F0RMAT(F10 . 7,5X ,E17. 10) 

85  FORMAT (F10. 7, 5X,’  NOT  COMPARED’) 

C 

RETURN 

END 
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